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ABSTRACT 


This  thesis  develops  methods  of  frequency  analysis  and  synthesis 
of  digital  computer  programs  describable  in  the  form  of  a linear  difference 
equation  with  constant  coefficients. 

The  mainspring  of  this  investigation  was  the  need  for  dealing 
with  control  systems  consisting  of  both  analog  and  digital  filters. 

Most  conventional  control  systems  consist  Of  analog  units  and  operate 
on  continuous  data,  but  digital  computers  use  sampled  data.  A uniform 
treatment  of  the  two  types  of  data  is  essential  in  the  analysis  of  control 
systems  incorporating  a digital  computer.  The  conventional  method  of 
treating  systems  operating  on  only  continuous  data  uses  Fourier  or  Laplace 
transformation!  that  is,  transformation  to  the  frequency  domain.  The 
conventional  method  of  treating  digital  programs  is  numerical  analysis,  which 
deals  almost  exclusively  in  the  domain  of  the  independent  variable!  that 
is,  the  time  domain.  By  exploiting  and  further  developing  those  areas  of 
numerical  analysis  to  which  frequency-transformation  techniques  were 
applied,  the  thesis  points  the  way  to  a common  language  of  dealing  with  a 
mixed-data  system. 

If  dataware  sanpled.  at  equal  intervals  of  time  (a  practical 
feature),  description  of  a linear  computer  program  always  reduces  to  a 
difference  equation.  It  is  possible  to  describe  such  a pregram  by  a transfer 
function  in  the  frequency  domain  in  a manner  analogous  to  the  conventional 
description  of  analog  filters.  Whereas  components  using  continuous  data 
have  transfer  functions  which  are  rational  functions  of  the  complex  frequency 
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variable  s,  those  of  a digital  program  are  rational  functions  of  % ■ e , 
where  e is  the  Naperian  base  and  T is  the  constant  interval  of  sampling. 

Having  described  the  digital  computer  with  its  program  by  a 
transfer  function,  one  may  apply  all  the  techniques  of  complex-variable 
and  transform  theory  to  deal  with  digital  filters.  Theorems  on  realizability, 
stability  and  other  properties  of  programs  are  developed,  and  the  amplitude, 
phase  and  locus  of  a program  are  defined.  The  adaptation  of  the  methods 
of  analog  filters  to  digital  ones  is  direct,  although  the  necessary 
modifications  are  often  significant. 

The  synthesis  of  computer  programs  can  be  conducted  along  lines 
employed  in  the  synthesis  of  networks.  First,  the  desired  frequency  charac- 
teristics of  the  program  are  stated;  next,  a rational  function  of 
-sT 

z ■ e is  found  which  approximates  the  desired  characteristics  for  real 
frequencies,  s » J«;  finally,  the  program  is  realized  on  basis  of  the 
approximating  transfer  function.  For  facilitating  the  approximation 
basic  entities  or  blocks  of  programs  are  analysed  and  methods  are  shown 
by  which  such  progranning  units  can  be  combined  to  obtain  the  frequency 
characteristics  of  the  complete  program.  Various  methods  of  program 
realization,  that  is,  programming,  are  developed  and  compared  on  the 
basis  of  time  and  storage  requirements,  and  criteria  are  developed  to 
permit  the  choice  of  the  optimum  programming  procedure  by  considering 
the  mere  form  of  the  program  transfer  function. 

Numerous  examples  of  program  analysis  and  synthesis  are  shown, 
and  one  example  of  synthesizing  a progran  for  the  compensation  of  a control 
system  is  worked  out.  The  latter  example  shows  that  the  frequency  analysis 
of  a complete  hybrid  system  can  be  undertaken  along  the  conventional  lines 
and  that  digital  compensation  of  a control  system  is  possible. 

ir 
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The  application  of  the  methods  of  the  thesis  to  various  problems 
in  numerical  analysis  is  also  shown.  The  problems  of  convergence  (stability) 
and  of  truncation  errors  (approximation)  can  be  analysed  in  the  frequency 
domain  effectively.  The  study  of  convergence  $y  conformal  mapping  is  related 
to  the  usual  methods,  and  a novel  way  of  estimating  truncation  error  is 
shown  provided  only  that  the  function  to  which  the  numerical  process  is 
applied  can  be  described  by  its  frequency  spectrum. 
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INTRODUCTION 

The  use  of  digital  computers  in  control  systems  is  now  coming 
into  the  fore.  Unlike  most  conventional  control  systems  involving  analog 
units  which  operate  on  continuous  data,  a control  system  employing  a 
digital  computer  of  the  present-day  type  must  use  sampled  data  in  the 
part  of  the  system  involving  the  digital  computer.  Hence,  some  parts 
of  this  system  use  continuous  data  and  others,  sampled  data.  The  Fourier 
and  Laplace  transform  methods  of'  analysing  continuous -data  control  systems 
is  well-known  and  developed,  but  the  conventional  treatment  of  digital 
computer  programs  is  by  numerical  analysis  or  in  the  time  domain.  There- 
fore, in  order  to  apply  the  methods  of  frequency  analysis  to  control 
systems  involving  digital  computers  (mixed-data  systems),  the  sampled 
data  part  of  the  system  must  be  described  in  the  frequency  domain.  Some 
work  along  these  lines  has  been  done  but  it  must  be  further  developed. 

An  analog  system  is  a physical  model  of  a set  of  differential 
equations j whereas,  a digital  system  is  a physical  model  of  a set  of 
difference  equations.  Operational  and  transform  methods  have  been  applied 
to  difference  equations  for  somB  time.  In  19l£  Gardner  and  Barnes1 
presented  a comprehensive  and  systematic  treatment  of  the  solution  of 
linear  difference  equations  with  constant  coefficients  by  the  Laplace 
transform  method.  However,  they  do  not  deal  with  stability  and  errors 


T 


Gardner  and  Barnes,  Transients  in  Linear  Systems,  John  Wiley  and  Sons 
New  Tork,  19  U2,  Chapter  It. 
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which  are  important  in  control  applications.  The  control  point  of  view 

is  stressed  in  Tustin's^  work  on  time  sequences.  In  19U9  and  1950  TustLn's 

2 

method  was  further  developed  by  Madwed  , who  shows  the  relations  of  his 
aspects  of  stability,  but  they  do  not  analyze  the  errors  associated  with 
their  approximations. 

3 

In  the  meantime,  Hurewicz  pioneered  the  analysis  of  pulsed 
filters  in  the  frequency  domain,  developed  stability  criteria,  and 
showed  several  examples  of  choosing  parameters.  It  should  be  noted, 
however,  that  Hurewicz* s filters  are  only  simple  units  such  as  differ- 
entiators and  lead  networks,  which  are  incapable  of  performing  involved 
computations  as  a computer  can.  Also,  Hurewicz  evaluates  the  output  of  a 
pulsed  filter  at  the  sampling  instants  only.  The  behavior  of  the  filter 
between  pulses  remains  a separate  problem,  and  no  ready  method  is  pre- 
sented to  investigate  the  whole  question  in  the  frequency  domain. 

W.  K.  Linvill^  shows  that  sampling  a continuous  function  is 
equivalent  to  the  modulation  of  a series  of  unit  impulses  by  the  function. 
The  result  is  a new  time  function  which  can  be  thought  of  as  being  applied 
to  the  sampled  data  part  of  the  system.  Furthermore,  this  new  time 
function  has  a Laplace  transform*  thus  a frequency-domain  analysis  is 
possible.  Linvill  shows  that  reconversion  from  discontinuous  to  continuous 


Tustin,  A Method  of  Analysing  the  Behavior  of  Linear  Systems  in 
Terns  of  Tias  Series,  J.I.EfcT  Vol.  9U,  Part  2a,#1,  pp.  130  - Ih2. 

2 

Madwed,  Number  Series  Method  of  Solving  Linear  and  Non-Linear 
Differential  Equations,  SC.D.  thesis  in  Mechanical  Engineering,  MIT. 

Hurewicz,  Filters  and  Servo  Systems  with  Pulsed  Data,  Chapter  5 of 
James,  Nichols  and  Phillips,  Theory  of  Servomechanisms. 

^ Linvill,  W.K.,  Analysis  and  Design  of  Sampled-Data  Control  Systems, 
Digital  Computer  Laboratory,  MIT,  Report  R-170. 
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data  is  a filtering  process  and  also  shows  what  happens  when  the  loop 
is  closed  on  a mixed  data  system.  He  is  concerned  only  with  the  effect 
of  sampling  on  the  system  and  does  not  consider  the  influence  of  digital 
computer  operations  on  the  system. 

This  report  is  a summary  of  the  work  done  by  Salzer^.  His 
results  permit  the  analysis  of  linear  digital  computer  programs  in  the 
frequency  domain j i.e.,  the  operation  of  a digital  computer  program  is 
described  by  a transfer  function.  Thus  the  field,  is  opened  for  the 
conplete  analysis  and  synthesis,  wholly  in  the  frequency  domain,  of  control 
systems  employing  digital  computers. 

From  the  frequency-domain  point  of  view,  conditions  governing 
the  realizability  of  program  transfer  functions  are  developed,  the  problem 
of  stability  is  studied,  and  conditions  to  insure  stability  are  given. 

Three  methods  of  realization  of  programs  from  their  transfer  functions  are 
presented,  and  the  time  and  storage  requirements  of  each  are  studied.  An 
elementary  example  of  transfer  function  synthesis  is  given.  As  in  the  case 
of  network,  theory,  the  analysis  of  a computer  program  in  the  frequency  domain 
is  straightforward  with  a -unique  result,  but  the  synthesis  of  a transfer 
function  has  many  alternate  realizations.  Also  as  in  network  theory,  the 
characteristics  of  the  transfer  function  to  be  realized  may  not  be  given 
directly  in  a forai  leading  to  immediate  realization  but  an  intermediate 
approximation  problem  may  need  to  be  solved.  Die  background  for  solving 
the  approximation  problem  has  been  set  up  in  that  conditions  of  physical 

Salzer,  J.M.,  Treatment  of  Digital  Control  Systems  and  Numerical  Processes 

in  the  Frequency  Domain,  SC.D.  thesis  in  Electrical  Engineering 
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realizability  have  been  derived  and  methods  of  realization  of  all 
realizable  transfer  functions  have  been>  obtained.  While  some  work  has 
been  done  directly  on  the  approximation  problem,  much  remains  to  be  done 
in  this  respect. 

The  function  of  this  report  is  to  provide  a concise  picture 
of  the  frequency  analysis  of  digital  control  systems  and  numerical  pro- 
cesses. The  first  chapter  describes  the  processes  of  sampling  and  de- 
sampling  continuous  functions  and  indicates  that  sampling  is  analogous 
to  impulse  modulation  while  desampling  is  analagous  to  ripple  filtering 
in  demodulation.  Thinking  of  sampling  as  impulse  modulation  allows  one 
to  relate  the  sanpled  to  the  continuous  function  in  either  the  frequency 
domain  or  the  time  domain.  Furthermore,  thinking  of  sampled  functions 
as  impulse  modulated  functions  allows  one  to  characterize  linear  computer 
operations  on  the  sampled  functions  by  transfer  functions. 

Chapter  II  derives  the  conditions  of  physical  realizability 
for  computer-program  transfer  functions,  discusses  stability  conditions 
on  these  transfer  functions,  and  presents  procedures  for  plotting  transfer 
loci. 

Chapter  III  deals  with  techniques  for  realization  of  transfer 
functions  with  some  attention  to  the  approximation  problem,  while  Chapter 
IV  deals  with  frequency  analysis  of  some  numerical  integration  formulas. 
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CHAPTER  I 

DESCRIPTION  OF  THE  SAMPLING  PROCESS 

1.1  Sampling  a Continuous  Function 

A digital  computer  operates  on  numbers  that  represent  samples  of 
continuous  signals  taken  at  discrete  instants  of  time*  The  time  interval, 

T,  between  samples  is  a constant  as  shown  in  Figure  1*1,  page  7 . In  this 
case,  the  input  to  the  computer  is  the  sampled  function,  i (t).  Hie  com- 
puter senses  the  amplitude  of  each  of  these  pulses  (as  a number)  and 
operates  on  the  number. 

Hie  purpose  of  this  chapter  is  to  describe  the  sampling  process, 
to  characterize  it  mathematically,  to  evaluate  how  well  a continuous  signal 
may  be  represented  by  its  saitples,  and  to  show  how  and  under  what  conditions 
a continuous  signal  may  be  recovered  from  its  samples. 

Hie  mathematical  model  of  the  sampling  process  which  will  be  de- 
rived later  is  very  similar  to  actual  physical  processes.  For  example, 
assume  that  i (t)  is  the  voltage  across  a pair  of  terminals  of  some  net- 
work. How  might  it  be  sampled?  The  voltage  may  be  sampled  by  connecting 
a condenser  across  the  terminals,  allowing  a current  flow  to  build  up  a 
charge  on  the  condenser  until  the  condenser  voltage  is  equal  to  the  terminal 
voltage,  and  then  disconnecting  the  condenser.  In  order  that  the  condenser 
voltage  be  equal  to  the  terminal  voltage  at  some  instant  of  time,  the 
sampling  time  should  be  as  small  as  possible.  It  may  be  made  very  small,  but 
not  zero.  Hie  total  charge  on  the  condenser  is  the  integral  of  the  current 
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flowing  over  the  time  required  to  take  the  sample.  Thus,  as  the 

sampling  time  decreases,  the  current  intensity  must  increase. 

Physically  this  is  how  sampling  might  be  done.  Ideally,  however, 
we  wish  to  take  the  sample  instantaneously  or  in  zero  time.  Therefore, 
for  ideal  sampling  in  the  above  example  the  current  flow  must  be  infinite 
for  zero  time  at  each  sampling  instant.  Thus,  in  the  ideal  case  the 
charging  current  is  an  impulse  whose  area  equals  the  amount  of  charge 
required  to  build  up  the  condenser  to  the  sampled  value.  Physically, 
ideal  sampling  is  not  possible,  but  the  idea  permits  us  to  set  up  a 

i 

model  of  sampling  that  can  be  treated  mathematically. 

1.2  Equivalent  Mathematical  Model  of  Ideal  Sampling  - Impulse  Modulation 

The  ideal  situation  in  the  above  example  is  to  transfer  to  the 

plates  of  the  condenser  a portion  of  charge  in  zero  time,  or  to  "hit"  the 

condenser  with  an  impulse  of  current.  The  same  end  can  be  obtained  if  we 

modulate  the  voltage  waveform  with  an  infinite  series  of  unit  impulses 

separated  by  equal  intervals,  T,  as  shown  in  Figure  2.  The  area  of  any 

one  of  the  modulated  impulses  equals  the  value  of  the  input  function  at  the 

corresponding  instant  of  time.  Thus,  impulse  modulation  is  analogous  to 

1 

the  process  of  sampling.  The  samples  of  Figure  1.1  have  finite  height,  zero 
width,  and  zero  area;  therefore,  -the  sampled  function  does  not  have  a Laplace 
transform.  The  impulses  of  Figure  1.2^  have  infinite  height,  zero  width, 


The  bar  (-)  over  i(t)  indicates  the  sampled  functions. 

The  circumflex  over  i(t)  indicates  the  impulse  - modulated  function. 
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but  non-zero  area;  therefore,  the  impulse -modulated  function  does  have  a 
Laplace  transform,  which  is  why  this  mathematical  model  has  been  set  up. 

1.3  Use  of  Impulse  Modulated  Functions  in  the  Analysis  of  Linear  Digital 
Computer  programs 

A digital  oomputer  operates  on  numbers  that  occur  at  discrete 
instants  of  time,  i.e.  it  operates  on  samples  of  a continuous  function. 

In  the  previous  section  it  was  rfiown  that  for  the  ideal  case,  sampling  is 
equivalent  to  impulse  modulation.  If  we  think  of  the  computer  as  "sensing” 
the  amplitude  of  samples,  we  may  just  as  easily  think  of  it  as  "sensing" 
the  area  of  impulses.  With  this  extension  or  mathematical  model,  we  may 
analyze  computer  programs  by  describing  the  input  to  the  computer  as 
impulses  instead  of  samples.  Since  a sample  does  not  have  a Laplace  trans- 
form, while  an  impulse  does,  the  advantage  of  this  extension  is  immediately 
obvious.  In  this  mathematical  model,  both  input  and  output  are  treated 
as  impulses,  and  both  have  Laplace  transforms.  In  conventional  (continuous- 
data)  systems,  the  transfer  function  is  the  ratio  of  the  transform  of  the 
output  to  that  of  the  input.  Since  both  input  and  output  of  computer 
programs  (when  treated  as  impulse-modulated  functions)  have  transforms, 
we  may  define  the  transfer  function  of  a linear  computer  program  as  the 
ratio  of  the  transform  of  its  output  to  the  transform  of  its  input.  In  order 
to  carry  out  this  analysis,  we  must  have  a knowledge  of  some  of  the  properties 
of  impulse -modulated  functions,  or  impulsed  functions.  The  remainder 
of  this  chapter  is  devoted  to  a discussion  of  some  of  these  more  useful 
properties. 
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l.h  Laplace  Transforms  of  Impulse-Modulated  Functions^- 

OUr  analysis  of  computer  programs  is  restricted  to  the  cases  in 
which  the  time  interval  between  samples  is  a constant,  T.  Thus,  the 
impulsed  function  can  be  expressed  as  the  product  of  a continuous  input 
function  and  an  infinite  string  of  unit  impulses,  the  interval  between 
impulses  being  T. 

As  the  following  derivation  will  show,  the  process  of  impulse 
modulation  may  be  readily  described  in  the  frequency  domain.  Essentially, 
since  the  string  of  unit  impulses  (which  is  the  carrier)  has  all  harmonics 
of  equal  amplitude,  the  impulse  modulated  wave  has  an  infinite  number  of 
side-bands  rather  than  just  the  two  which  are  present  for  a sinusoidal 
carrier.  The  method  of  the  derivation  is  to’make  a ..Fourier'  analysis'  of 

the  carrier  and  to  associate  each  side-band  of  the  impulse  modulated  wave 
with  a Fourier  component  of  the  carrier.  Let  i(t)  be  the  continuous  input 
function  and  “ kT)  be  the  infinite  string  of  unit  impulses. 

(x)  = unit  impulse  occurring  at  x = 0.]  Then  the  impulse -modula te d input 
function  is, . 

i(t)  = i(t)  ^-=-  kT).  (1-1) 


To  find  the  Laplace  transform  of  (1-1)  let  us  first  find  the 
complex  Fourier  series  of  the  string  of  unit  inpulses. 


H(t 


k = 


- kT) 


m = - 1 


c e 
m 


jmjfLb 


(1-2) 


4" 


A more  couple  te  deriva  tion  and  discussion  of  the  transforms  of  impulse 
modulated  functions  is  given  in  Reference  2,  page  3. 
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In 


d-2), -n-  = V-  . 


Ihe  c _ s are  the  complex  Fourier  coefficients, 
m 


Solving  for  c in  the  usual  manner  we  have. 

IQ 

T/2 


ST  “I 

m * ? 1 f ^ " kT>  dt.  (1-3) 

J L f^TT J 


-T/2 


Qjr  writing  out  a few  terms  of  the  series.  (1-3)  becomes.. 


T/2 

>B  = | j I + |l(t  - T)  + n(t)  + n(t  + T)  l...j6'>Xltdt 

-T/2  L 


Cl-lt) 


Within  the  range  of  the  integral,  the  only  term  inside  the  bracket  of  the 

integrand  that  is  non-zero  is  the  term,  p,( t).  Ulus  (l-li)  becomes. 

T/2 


5«s!  J 

-T/2 


n(.t)  e“^m  ^ dt. 


(1-5) 


Because  of  the  unit  impulse  in  the  integrand,  the  value  of  the  integral  is 
just  evaluated  at  t = 0,  which  is  unity.  Therefore, 


cm  ~ 


(1-6) 


and  the  Fourier  series  of  a string  of  unit  impulses  is, 

im-TLt 


^(t-JcT)  * | 


k = -* 


(1-7) 


a = -« 
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Then  the  impulsed  function  becomes, 


T(t>-  4a 


Jm-fLt 


m = - — ^ 


(1-8) 


Now  takB  the  Laplace  transform  of  the  above  equation. 


'l(i)  = L JjL(t) 

1- 

i?~~~  ejmiU  1 

m = -^o  -* 

• (1-9) 

The  indicated  summation 

can  be  done 

after  the  transformation 

is  made. 

~(s)  = | 

0*7 

> 

m = - 

L j^i(t)  e^^j 

(1-10) 

A fundamental  theorem  in  Laplace  transform  theory  leads  directly  to  the 
following  result: 

^(=)=i 

(s  + jmXL) 

- (1-11) 

^ 1 

m = - 


Thus  we  see  that  the  Laplace  transform  of  an  impulsed  function  is  periodic 

having  a repetition  interval  of  j-O-. 

/\/ 

An  important  fact  about  I (s)  should  be  observed  from  (1-11).  It 
is  that  there  is  a unique  correspondence  between  I (s)  and  I (s)  if  and 
only  if  the  frequency  spectrum  of  i(t) } the  continuous  time  function,  lies 
in  the  Irange,  If  the  spectrum  of  i(t)  lies  outside  this 

range,  ^(s)  will  specify  the  spectrum  (in  the  range  <t>  of  a 

continuous  time  function,  but  this  time  function  will  differ  from  the 


I 
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f 

a 


original  time  function.  Thus,  there  is  a limitation  of  bandwidth 
caused  by  sampling.  Figure  1,3  illustrates  the  case  of  unique  corres- 
pondence, and  Figure  i.i*,  the  case  in  vhich  the  spectrum  of  i(t)  is  too 
wide. 


As  given  by  (1-11),  l(s)  consists  of  an  infinite  number  of  terms} 
however,  an  infinite  series  is  difficult  to  handle,  and  it  is  desirable 
to  have  a closed  form  expression  for  I(s).  This  can  be  obtained  from 

K. 

the  partial  fraction  expansion  of  I(s).  Consider  a -typical  term, 


s - s^* 


of  the  partial  fraction  expansion  of  I(s).  Referrihg  to  (1-11)  we  see 
that  corresponding  to  this  typical  term,  'T(s)  will  have  a typical  series  of 
terms  of  the  form. 


1 

> ~s  - s1  + jk  * 


Thus  we  see  that  the  pole  at  s = s.  is  repeated  ap.  infinite  number  of  times 
at  intervals  of  j the  line  through  these  poles  being  parallel  to  the 


imaginary  axis  in  the  s - plane. 


A closed  form  equivalent  of  the  above  typical  series  can  be 
obtained  by  a change  of  variable  in  the  following  equation.1 


tt  fc-  cot  ir  i = 1 + 2 


HI—  -g  _ 2 (1-12) 

n = 1 n 


1 

Knopp,  "Theory  and  Application  of  Infinite  Series",  New  York,  19  U8,  p.  1*19 


C.  Spectrum  of  ij(t)  that  would  produce  the  same  sampled  function 
as  (B). 


Figure  l.U  Illustration  of  Bandwidth  Limitation  Caused 
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Divide  each  side  of  (1-12)  by  and  make  the  change  of  variable,  & » j — . 


»“*«  5r  = tr*  jt  1 


oa 


Go 


n = 1 


JCT  2 (W3) 

E*  + n 


Multiply  both  sides  of  (1-13)  by  j/fu  and  obtain, 


. n . . rr  as  _ n ..  tt  oo  _ 1 

cot  j -fj-  - JJ-  coth  3C-  - s ♦ 5 


2co 


ry  os2  + n2_TL  2 


(1-lU) 


The  infinite  series  of  poles  of^I(s)  corresponding  to  a pole  of 
I(s)  at  s.^  can  be  put  into  a form  that  is  identical  to  the  right-hand  member 
of  (l— 1I4.)  as  follows:  separate  the  term  for  k = 0. 


s - si  ♦ jk  JL  s 


” si  + ^Fa"l^"Si  + s~si  “ 


Combine  -the  two  terms  in  the  summation. 

1 - 1 . w 

s - si  + jkJl.  s - y 

^ k = k = 1 


2(g  - 3p 

(s  - si)2  + k2  -Tt2 


(1-15) 
a-16) ) 


1 comparison  of  (l-ll:)  and  (1-16)  shows  that. 


■£= 


_i = * cott  *■  (■  - H). 

- s±  + jkA.  THT  _TU 


(1-17) 


Thus  we  have  the  following  closed  form  equivalent  of  the  typical  series  of 

r(a), 

1 t Kj  m 

2. »-».«!  JlcJL  = T coth  5 (s  - .j), 

ks.oo  1 


Ki 

T 


(1-18) 
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fbr  a pole  of  I(s)  at  s = s^.  Therefore,  corresponding  to  the  partial 
fraction  expansions  of  l(s),  we  have  the  following  series  for^s). 

n 

“iCs)  = ^ y ~~  K±  coth  | (s  - s^  , 
i = 1 

where  "n"  is  the  total  number  of  poles  of  l(s),  taking  into  account  multiple 
«> 

poles. 


Let  us  now  investigate  the  limitations  on  the  positions  of  the 
poles  of  I(s)  due  to  sampling.  Consider  an  infinite  strip  of  width  -TV.  in 
the  a-plane  and  parallel  to  the  real  axis  as  shown  in  Figure  1.5.  Assume 
that  aLl  the  poles  of  1(a)  lie  within  this  strip  and  in  the  left  half  plane 
(LHP).  Thus,'l(s)  has 


jw 


Figure  1.5  Infinite  Strip  Containing  Poles  of  I(s) 
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poles  at  these  points  plus  poles  at  points  shifted  from  the  s.^  ' s and 

s*^  's  (#  means  conjugate)  by  the  distance  - jk  _TL  . Sinoe  l(s)  has 

poles  only  in  the  strips  being  considered,  there  is  a one-to-one  correspondence 

between  the  poles  of  I(s)  andTf(s)  that  lie  in  the  same  strip.  However, 

if  I(s)  had  poles  outside  this  strip,  there  no  longer  would  be  this  one-to-one 

correspondence. 
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CHAPTER  II 

TRANSFER  FUNCTION  OF  COMPUTER  PROGRAMS  - REALIZABILITY  AND  STABILITY 


Using  the  properties  of  impulse -modulated  functions  given  in 
Chapter  I,  we  are  now  ready  to  investigate  transfer  functions  of  computer 
programs.  Our  interest  in  program  transfer  functions  is  much  more  than 
academic.  The  transfer  function  describes  the  program  completely  and 
with  it  we  car,  analyze  and  synthesize  control  systems  employing  digital 
computers  by  conventional  frequency  domain  methods. 

In  this  chapter  a linear  digital  computer  program  is  defined  in 
terms  of  the  mathematical  model  of  sampling  set  up  in  Chapter  I,  its  transfer 
function  is  derived,  and  methods  for  dete mining  the  realizability  and 
stability  of  transfer  functions  are  given.  Several  examples  of  stability 
determination  are  also  presented. 

2.1  Transfer  Function  of  Linear,  Real-Time,  Digital  Computer  Program 
As  pointed  out  in  Chapter  I,  the  input  to  a digital  computer 
may  be  assumed  to  be  an  impulsed  function,  for  purposes  of  mathematical 
analysis.  A linear  program  of  a digital  computer  operating  in  real  time 
is  one  in  which  the  present  output  is  a linear  function  of  the  present 
and  past  inputs  and  the  past  outputs.  The  general  form  of  this  relation 
is, 

~o(t)  -y~~  s^TCt  - kT)  - 5"°—  bk  <3r(t  - kT),  (2-1) 
k • 0 k - 1 


in  which  all  a^'s  and  b^'s  are  real,  and  T is  the  time  between  samples « 
The  time  required  for  the  computation  must  be  less  than  T if  each 
calculation  is  to  be  completed  before  the  next  input  arrives. 
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Taking  the  Laplace  transform  of  (2-1)  yields, 

oc.)  - 1(.)  y-  s .-k*T yooy8-  bk  .-*»*.  o>-2) 

k - 0 k - 1 

As  in  continuous  data  systems,  we  will  define  the  transfer  function 
of  a computer  program  as  the  ratio  of  the  transform  of  the  output 
to  that  of  the  input.  Let  W(s)  be  the  transfer  function  of  a computer 
program}  then, 

W(s)  - • (2-3) 

1(e) 

Solving  for  T5ts)/  its)  from  (2-2)  we  obtain. 


W(8)  * 


fr(») 

T(a) 


(2-U) 


as  the  transfer  function  of  a linear,  real-time,  digital  computer  program. 
With  the  understanding  that  bQ  ■ 1,  (2-U)  becomes, 

m 


W(s) 


ZZ\ 


-ksT 


-ksT 

e 


(2-5) 


k » 0 


The  inverse  steps  from  (2-5)  to  (2-l)  are  unique;  therefore, 

(2-5)  is  the  general  form  of  the  transfer  function  of  a realizable,  linear, 
digital  computer  program, . Thus,  to  be  realizable,  the  transfer  function 
of  a linear,  digital  computer  program  must  be  expressible  as  the  ratio  of 
two  polynomials  in  e • The  criteria  for  stability  will  be  discussed  in 


a later  section 


Report  R-225 


-20. 


It  has  already  been  shown  that  the  Laplace  transform,  l(s),  of 
the  impulsed  input  function  is  periodic  of  period as  seen  in  (l-ll). 

By  showing  that  W(s)  is  also  periodic  with  the  same  period,  we  can  prove 
that  ots)  is  also  periodic  of  period  -A. . A typical  term  of  either  numerator 

i._m 

or  denominator  of  W(s)  contains  e“  . For  s — >s  ♦ jm-fl  (m  is  a positive  or 

negative  integer),  the  typical  term  becomes, 

-k(s  ♦ jm/L)T  -ksT  -jkm-flT 
6 * 6 8 • 

% 

As  T-O.  - 2ir  and  k and  ra  are  integers,  the  second  factor  is, 

-jkra  /U  -j2fikm  , 

S ■ 0 • X • 

u -k(s  ♦ jm/l)T  -ksT 

Therefore,  the  terms  of  the  numerator  and  denominator  Df'  W(s)  are  periodic 
of  period/"!,  and  so  is  W(s).  In  equation  form  this  means,  W(s)  » 

¥(s  ♦ jm/1),  for  m a positive  or  negative  integer.  The  product  of  two 
periodic  functions  is  also  periodic.  Since  0(s)  - W(s)Tl(s),  ^(s)  is  also 
periodic  of  period_fl,  as  indeed  it  should  because  the  computer  output  is 
also  sampled. 

Since  all  the  coefficients  of  (2-5)  are  real,  it  is  readily 
* * 

seen  that  Vf(s)  - W(s  ) , in  which  the  asterisk  means  conjugate.  For  real 
frequencies  this  becomes  W(jeo)  - W(-jco)  . This  fact  together  with'  the 
periodicity  of  W(s)  tells  us  that  W(s)  is  completely  specified  for  all  s 
if  it  is  defined  over  the  range,  0 

Suramaryi  In  order  to  be  realizable,  the  transfer  function  of 
a linear  digital  computer  program  must  be  expressible  as  the  ratio  of 

-Ip 

polynomials  in  e • W(s)  is  periodic  of  perLod^fl } io6»#  W(s)  ■ 

V(s  ♦ jm/1).  Specification  of  W(s)  over  the  range,  0 4=  oa  completely 
determines  V(s). 
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2.2  Stability  of  Programs 

We  have  expressed  the  transfer  function  of  computer  programs  as 
a function  of  the  complex  frequency  "s"|  therefore,  the  same  methods  of 
investigating  stability  as  used  in  network  analysis  and  servomechanisms 
are  applicable.  The  general  necessary  and  sufficient  criterion  for  stability 
of  a unit  is  that  its  transfer  function  have  no  poles  in  the  right  half 
s-plane  (RHP)  or  multiple  poles  on  the  jonaxis.  In  network  analysis  the 
frequency-domain  method  used  to  study  stability  is  to  map  a contour 
enclosing  the  right  half  of  the  s-plane  (the  contour  is  usually  the  jco-axis 
and  an  infinite  semicircle)  into  the  W-plane.  Because  of  the  transcendental 
nature  of  the  transfer  function  of  a realizable  computer  program,  the 
mapping  contour  in  the  s-plane  must  be  modified. 

r 

As  we  have  shown  before,  the  transfer  function  of  the  computer 


program  is. 


V(s)  “ PCs)  - 

QCs) 


ZZ> 

k - 0 
n 

k - 0 


-ksT 

e 


-ksT 

e 


(2-6) 


in  which  P(s)  is  the  numerator  and  Q(s),  the  denominator  of  W(s);  and  it  is 
assumed  that  P(s)  and  Q(s)  have  no  common  factor.  Both  P(s)  and  Q(s)  are 
entire  transcendental  functions  having  as  their  only  singularity  an  essential 
singularity  at  infinity. t Hence,  we  see  that  the  only  singularities  of 
W(s)  in  the  finite  s-plane  are  poles,  and  these  poles  occur  at  the  zeros 
of  Q(s).  Our  stability  criterion  is  that  there  be  no  poles  of  W(s)  in  the 
RHP  and  only  simple  poles  on  the  imaginary  axis.  Therefore,  in  order  for 


t For  a further  discussion  of  entire  transcendental  functions,  consult  Knopp, 
"Theory  of  Functions,"  or  Guillemin,  "The  Mathematics  of  Circuit  Analysis." 
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the  pro gran  to  be  stable,  Q(s)  must  have  no  zeros  in  the  RHP  and  only  staple 
zeros  on  the  joa-axLs . To  investigate  the  possibility  of  Q(a)  having  zeros 

in  the  RHP  or  on  the  imaginary  axis,  we  may  take  advantage  of  the  periodicity 
of  Q(s).  In  proving  that  W(s)  is  periodic,  it  was  shown  that  e"  is 
periodic  property.  Therefore,  if  Q(s)  has  a zero  in  the  RHP,  it  must  have 
one  in  the  send-inifinite  strip  shown  in  Fig.  2d. 


Figure  2.1  Semi -ini finite  strip  of  s-plane  that  must  have  a zero 
of  Q(s)  if  Q(s)  has  any  zeros  in  the  RHP, 

-sT 

Consider  the  map  of  the  contour  of  Fig.  2.1  into  the  e plane. 
Let  us  begin  the  path  at  the  origin  in  the  s-plane  and  encircle  the  strip 
in  a clockwise  direction,  corresponding  to  increasing  frequency.  It  is 

-sT 

readily  understood  that  corresponding  path  and  enclosed  region  in  the  e 
plane  is  shown  in  Fig.  2.2.  The  origin  of  the  s-plane  maps  into  the 
point  (1,0)  in  the  e”  plane.  The  corresponding  sections  of  the  path  are 
marked  by  small  letters  on  both  contours.  In  Fig.  2.2  we  see  that  the 
paths  (b)  and  (d)  cancel  leaving  the  annular  ring  as  the  region  conformal 
to  the  strip  of  the  s-plane  that  is  under  consideration.  As  0^  (of 
Fig.  2.1)  approaches  oo , the  radius  of  the  circular  path  (c)  in  Fig.  2.2 
approaches  zero.  Thus,  the  conformal  map  of  the  indicated  strip  consists 
of  two  separate  contours!  one,  a unit  circle  centered  at  the  origin 


Figure  2.2  Conformal  map  of  the  semi -inif ini te  strip  of  Fig.  2.1 
Into  the  e plane. 

and  the  other  an  infinitesimally  small  circle  that  excludes  the  origin 

in  this  particular  case.  Only  a slight  extension  of  the  foregoing  procedure 

_ST  .JrgT 

is  required  to  determine  the  map  of  powers  of  e . The  map  of  e 

■•sT 

will  appear  like  that  of  e except  that  each  of  the  two  separate  paths 

will  be  traversed  "k"  times/  the  region  excluded  by  the  infinitesimally 

small  circle  will  be  that  at  the  origin.  Thus  we  see  that  the  map  of 

this  semi -infinite  strip  of  the  s-plane  is  effective  in  handling  the 

essential  singularity  of  Q(s)  at  a©  . 

Now,  consider  the  conformal  map" of  the  semi -infinite  strip  of 

sT  — 2sT 

Fig.  2.1  into  the  Q-plane.  Remembering  that  Q(s)  * l*b,  e ♦ b^e  ♦ 

...  ♦ b^e"*118*,  we  see  that  the  map  of  this  strip  into  the  Q-plane  will 
exclude  the  point  (1,0),  (the  map  of  each  term  except  the  first  excludes 
the  origin).  This  eliminates  the  need  for  mapping  path  (c).  Moreover, 
since  the  paths  (b)  and  (d)  cancel,  we  need  to  plot  only  the  paths  (a) 
and  (e).  In  other  words  the  only  part  of  the  s-plane  contour  that  we 
need  to  plot  in  order  to  determine  the  locus  of  Q(s)  is  the  part  of  the 
contour  that  lies  on  the  imaginary  axis.  This  contour  in  the  Q-plane 
will  encircle  the  origin  z-N  times  in  the  counterclockwise  directipn,  where 
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s is  the  number  of  zeros  and  N is  the  number  of  poles  of  Q(s)  (taking  into 
account  their  multiplicity)  in  this  strip  of  the  s -plane.  It  has  already 
been  pointed  out  that  Q(s)  is  an  entire  transcendental  function  and,  therefore, 
has  no  poles  in  this  strip.  So,  N ■ 0,  and  the  contour  in  the  Q-plane  will 
encircle  the  origin -Z  times  (clockwise  is  to  be  understood).  The  condition 
for  stability  of  W(s)  states  that  Q(s)  must  not  have  any  zeros  in  the  RHP 
or  any  multiple  order  zeros  on  the  imaginary  axis.  Therefore,  the  map  in 
the  Q-plane  must  not  enclose  the  origin;  Z must  be  zero.  If  Q(s)  has  zeros 
on  the  imaginary  axis,  the  Q-plane  locus  will  pass  through  the  origin.  In 
this  case,  we  must  determine  the  order  of  the  zero.  The  following  method 
can  be  used*  Assume  that  the  locus  in  the  Q-plane  passes  through  the  origin 

for  s ■ Joo . . Then  Q(s)  must  contain  the  factor  (e  - e-J  i ) where 

1 2 
n is  the  order  of  the  zero.  Divide  Q(s)  by  (e”sT  - e”^6>iT)  . If  there  is 

no  remainder,  the  zero  is  of  higher  order  than  the  first  and  the  program 

will  be  unstable. 

In  addition  to  determining  the  stability  of  pro grans,  conformal 
maps  give  an  indication  of  the  degree  of  stability  or  instability  and  an 
approximate  value  of  the  frequency  at  which  the  program  is  or  may  become 
unstable.  The  amount  by  which  the  locus  in  the  Q-plane  misses  encircling 
the  origin  gives  a measure  of  the  stability  of  the  program.  The  farther  the 
locus  is  from  the  origin,  the  more  stable  or  convergent  the  program.  The 
frequency  corresponding  to  the  point  on  the  Q-plane  locus  nearest  the  origin 
is  approximately  the  frequency  at  which  the  program  is  or  may  become  unstable, 
or  at  which  it  will  oscillate  in  a damped  fashion.  , 


( ' > 


'I  ( 
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In  addition  to  expressing  the  program  transfer  function  as  a function 

-sT 

of  "s"  we  may  also  write  it  as  a function  of  e » Make  the  change  of  variable, 
-sT 

z - e o Then  we  may  define  a new  function, 

m i. 


> v 


(2-7) 


£ 


k - 0 


V 


V(z) 


aQ  ♦ a^z  ♦ a^z  ♦ 
_2 ? 

1 ♦ b^z  ♦ b^z  ♦ 


♦ a z 
m 


m 


(2-8) 


♦ b z 

n 


n 


It  is  readily  seen  that  the  right  half  of  the  s-plane  maps  into  the  inside 
of  a unit  circle  centered  at  the  origin  in  the  z-plane.  The  imaginary  axis 
of  the  s-plane  becomes  the  unit  circle  in  the  z-plane  (see  Fig.  2.3). 


Figure  2.3  Map  of  right  half  of  s-plane  into  z-plane 
Therefore,  if  the  program  is  to  be  stable,  all  the  zeros  of  D(z)  must  lie 
outside  the  unit  circle  except  that  single  order  zeros  may  occur  on  the 
unit  circle.  In  other  words,  the  magnitude  of  the  roots  of  D(z)  must  be 
greater  than  or  equal  to  unity,  and  the  roots  of  unity  magnitude  must  be 
simple. 
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Summary « To  test  for  the  stability  of  a program,  map  the  semi-infinite 

strip  of  Fig.  2.1  into  the  Q-plane  \Q(s)  is  the  denominator  of  the  program 

transfer  function^  If  the  locus  in  the  Q-plane  does  not  enclose  the  origin, 

the  program  is  stable  or  convergent.  If  the  locus  passes  through  the  origin, 

Q(s)  has  a zero  and  the  order  of  this  zero  must  be  determined.  If  the  zero 

is  of  first  order,  the  program  is  stable;  otherwise,  unstable.  An  alternate 

-eT 

method  isi  Make  the  change  of  variable,  z ■ e , and  find  the  magnitude 
of  the  z-roots.  If  each  root  has  either  a magnitude  greater  than  unity  or 
equal  to  unity  and  is  simple,  the  program  is  stable  or  convergent;  otherwise 
unstable  or  divergent. 

2.3  Loci  of  Q(s) 

In  the  previous  section  it  was  demonstrated  that  the  stability 
of  a program  can  be  determined  by  mapping  the  contour  enclosing  the 
semi -infinite  strip  of  Fig.  2.1  into  the  Q-plane.  It  was  also  shown  that 
the  only  part  of  this  contour  that  we  need  to  plot  is  that  on  the  imaginary 
axis.  The  paths  (b)  and  (d)  cancel  and  the  path  (c)  excludes  the  point 
(1,0)  in  the  Q-plane,.  Hence,  we  are  interested  in  the  properties  of  Q(ja>) 
and  its  locus  in  the  range, 

Q(ja>)  has  several  properties  that  are  helpful  in  determining  its 


10CU8 . 


(1)  Q(j<»)  is  periodic  of  period-O..  This  was  proved  in  the 
previous  section. 

(2)  Q(jo>)  * Q(-jco)  » This  property  follows  directly  from 
the  fact  that  Q(jco)  is  a polynomial  in  e”*^,  and  as  a 
consequence,  the  locus  of  Q(jco)  must  be  symmetrical  about 


the  real  axis 
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(3)  At  (n  - 0 and  * ^ Q(j<o)  is  real  and  its  locus  at  these  two 

points  crosses  the  real  axis  either  normally  or  tangentially. 
This  statement  is  proved  and  elaborated  upon  in  Appendix  A. 

As  a consequence  of  the  first  two  properties,  the  locus  of  Q(ja>) 
for  0 Jm  os  complet ely  determines  the  locus  in  the  Q-plane.  The  locus 

for  Z-0  is  Just  the  mirror  of.  that  for  the  positive  values  of  00. 

Thus  the  first  two  properties  result  in  a substantial  reduction  in  the  amount 
of  work  required  to  plot  the  locus  of  Q(jco).  The  third  property  enables 
one  to  determine  accurately  the  shape  of  the  locus  in  the  neighborhood 
of  00  ■ o and  . 

Several  methods  may  be  used  to  determine  the  locus  of  Q(jco). 

Three  of  these  arei  (l)  add  the  loci  of  the  individual  terms  of  the 
polynomial  (each  locus  is  a circle)  to  obtain  that  of  Q(ja>);  (2)  factor 

Q(jco)  and  multiply  the  loci  of  the  factors;  and  (3)  express  Q(ja>)  in  the 
form  R(<x>)/0(a>)  and  make  a point  by  point  plot.  The  method  that  is  best 
to  use  depends  on  the  particular  Q(jco).  However,  it  is  to  be  expected  that 
the  extra  analytic  work  required  in  methods  (2)  and  (3)  will  result  in 
less  graphical  work  and  more  accurate  loci.  None  of  the  methods  will  be 
discussed,  but  they  will  be  illustrated. 

Let  us  now  consider  several  examples  of  loci  of  QCjco). 

1.  Let  Q(s)  - 1 - 0.8  e'sT  ♦ 0.3  e"2sT  (2-9) 


Take  the  derivative  with  respect  to  s< 


- 0.8Te“sT  - 0.6Te"sT 
§§  - 0.2Te"sT  (U  - 3e’8T) 


(a) 


(b) 


(*-10) 


4 


1 ) 
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For  neither  s - 0 nor  — j^is  the  derivative  zero,  so  the  locus  is  normal 
to  the  real  axis  at  both  points.  Fig.  2.1*  (a)  shows  the  locus  of  Q(joo)  and 

O aw 

how  it  was  obtained  (for  the  point  at  - jj— ) from  the  loci  of  the  individual 
terms.  The  locus  for  negative  values  of  a is  shown  by  the  dashed  curve, 
since  it  may  be  obtained  from  the  other  half  of  the  locus.  After  this,  only 
the  locus  for  positive  a will  be  drawn.  The  locus  does  not  enclose  the  origin; 
therefore,  a program  whose  transfer  function  has  the  denominator, 

1 - 0.8e-8*  ♦ 0.3e-2s^,  will  be  stable. 


2.  Let  Q(s)  - 1 - 0.8e“sT  ♦ 0.1*e“2BT 


(2-11) 


-TV 


Take  the  derivative  with  respect  to  s. 

H - 0.8Te’sT  - 0.8Te"2sT 
- 0.8Te“sT  (1  - e“sT)  , 


(a) 

Cb) 


(2-12) 


For  s*0,  the  derivative  is  zero,  Q(s)  has  a saddle  point  here.  Q(s)  can 
be  rewritten  as, 

Q(s)  - 0.6  ♦ 0.1*  (1  - e*sT)  (2-13) 

which  brings  the  saddle  point  into  evidence.  In  this  case  p - 2,  so  at 
(o  * 0 . the  locus  is  tangent  to  the  real  axis . At  s - ♦ / 0,  so  the 

locus  is  normal  to  the  real  axis  at  this'  point.  Fig.  2.1*  (b)  shows  the 
resulting  locus.  It  does  not  enclose  or  pass  through  the  origin;  therefore, 
this  Q(s)  will  lead  to  a stable  progran. 

3.  Let  Q(s)  - 1 - 0.8e"sT  «■  0.$e"2sT  (2-UO 


Take  the  derivative  with  respect  to  s. 

g - 0.8Te”sT  - Te"28T 


- 0.2Te”aT  (1*  - 5e“sT) 


(a) 

(b) 


(2-15) 
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For  neither  8-0  nor  the  derivative  zero,  so  the  locus  of  Q(jco)  is 

perpendicular  to  the  real  axis  at  both  points.  The  locus  (as  shown  in 
Fig.  2*1}  (c ) ) does  not  enclose  the  origin;  therefore,  this  Q(s)  is  the 
denominator  of  a stable  program  transfer  function* 


Report  R-225 


Report  R-225 


CHAPTER  III 

ANALYSIS  AND  SYNTHESIS  OF  LINEAR,  DIGITAL  COMPUTER  PROGRAMS  IN  THE 

FREQUENCY  DOMAIN 

In  the  first  part  of  this  chapter  the  analysis  of  transfer  functions 
is  dealt  with  by  expanding  the  transfer  function  into  partial  fractions. 

Next,  programs  are  realized  from  transfer  functions  by  three  methods*  direct 
programming,  cascade  programming,  and  parallel,  programming;  and  the  storage 
and  time  requirements  of  each  are  presented*  In  the  last  part  of  the  chapter 
a short,  general  discussion  of  synthesis  is  given,  and  one  possible  synthesis 
procedure  is  illustrated  by  the  synthesis  of  a program  for  differentiation. 

3.1  Response  of  Programs  at  Real  Frequencies 

The  input  to  a computer  has  a certain  frequency  spectrum,  and 

in  order  to  analyse  the  action  of  a computer  program  on  this  input  function, 

we  need  to  have  a knowledge  of  the  frequency  response  of  the  program.  Thus, 

we  are  interested  in  the  locus  of  W(joo),  the  map  of  the  Jo-axis  of  the  s-plane 

into  the  W-plane.  A familiarity  with  the  frequency  characteristics  of  the 

simple  transfer  functions  is  essential  for  the  understanding  of  the  possibilities 

and  limitations  of  more  complicated  ones. 

In  many  cases  the  desired  locus  of  the  transfer  function  of  a 

digital  computer  program  is  given,  and  the  problem  is  to  approximate  this 

locus  by  that  of  a realizable  program;  i.e.,  by  a ratio  of  polynomials  in 
-T 

e • A study  of  the  loci  of  typical  terms  of  W(jo>)  is  helpful  in  making 
this  approximation. 

3.2  Analysis  of  Building  Blocks  of  Transfer  Functions 

As  we  have  seen,  the  transfer  function  of  a linear,  digital 
computer  program  is  most  generally  expressed  as  the  ratio  of  polynomials  in 
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WCs)  - Bisl  . a04  V * V 


-28T 


♦ a e 
m 


-kbT 


i(s) 


-i  . v -sT  -2sT  _ ^ -nsT 

1 ♦ b,e  ♦ b-e  ♦ ♦ a e 

id  n 


(3-1) 


A partial  fraction  expansion  of  W(s)  can  be  made,  and  we  may  caLl  the  individual 
terms  of  the  expansion  the  basic  building  blocks  of  a program’ transfer  function* 
In  general  W(s)  may  be  broken  up  into  a polynomial  plus  first  and  second  degree 
partial  fractions  (from  the  real  aid  conjugate  complex  roots  of  the  denominator, 
respectively) 

In  analysing  computer  programs  in  the  frequency  domain  [finding 
the  locus  of  W(jo#)J|  several  methods  can  be  used.  Two  of  these  arei  (l) 

Find  the  loci  of  the  numerator  and  denominator  polynomials  and  then  divide* 

(2)  expand  W(s)  into  partial  fractions,  find  the  locus  of  each  term  of  the 
expansion,  and  add  the  resultant  loci.  In  most  cases  the  first  method 
is  easier  to  use,  but  the  second  is  included  here  because  of  its  connection 
to  the  synthesis  of  program  transfer  functions  (approximation  of  a desired 
locus  by  a sum  of  the  basic  building  blocks).  A familiarity  with  some  of 
the  possible  loci  of  polynomials  and  first  and  second  degree  partial  fractions 
is  an  aid  in  the  synthesis  procedure. 

The  loci  of  second  degree  polynomials  have  already  been  discussed, 
and  the  locus  of  a fourth  degree  polynomial  will  be  illustrated  in  connection 
with  polynomial  building  blocks.  Since  there  is  only  a short  step  from 
the  loci  of  polynomials  to  the  locus  of  a transfer  function,  the  first 
method  of  finding  the  locus  of  W(jco)  will  not  be  discussed. 

For  use  of  the  second  method,  we  will  investigate  the  loci  of  typical 
terms  of  the  pertiel  fraction  expansion  of  W(j a)’. 
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3*21  Polynomials 

A typical  polynomial  transfer  function  is  of  the  form, 


For  s ■ jo>,  the  locus  of  each  term  of  the  polynomial  is  a circle*  For 

Q(s),  the  constant  term  is  unity,  but  for  polynomial  building  blocks,  the 

constant  term  may  have  any  real  value*  In  the  previous  chapter,  three 

-sT 

examples  of  the  locus  of  second  order  polynomials  in  e are  given,  so  nov 
let  us  find  the  locus  of  a fourth  order  polynomial.  Let, 

¥Cs)  - ^ (13  ♦ Ue"aT  - 3e"2aT  ♦ 2e"3sT  - e*4*8*)  (3-3) 

First  examine  the  function  for  saddle  points. 

U - g (-UTe”sT  ♦ 6T«'2sI  - 6To'3aT  ♦ liTe-**”1)  (3-ll) 

g - - g e-1  (2  - 3.-ST  * 3.'2,T  - 2.-3”1).  (3-5) 

The  derivative  is  zero  far  s - 0,  therefore  W(s)  has  a saddle  point  there* 
W(s)  can  be  written  in  the  form 

W(s)  - 1 - (2  ♦ e‘2sT)  (1  - e*sT)2  , (3-6) 

which  brings  the  saddle  point  at  s » 0 into  evidence.  The  saddle  point 
is  of  first  order;  therefore,  the  locus  is  tangent  to  the  real  axis  at 
s » 0.  W(jco)  is 

WCjco)  - 1 - (2  + e"J2wT)  (1  - e-J“T)  (3-7) 

In  this  case  it  is  easier  to  determine  the  locus  of  W(jaj)  by 
plotting  the  second  term  of  (3-7)  and  shifting  the  origin  one  unit  to  the 
left.  A convenient  way  to  find  the  locus  of  the  second  term  is  to  let 

- (2  ♦ e"J2aiT)  (1  - e"**1)2  - Re^  , (3-8) 
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If  the  constants  *4.  and  P are  real,  then  (3-11)  cai  be  considered 
a basic  building  block.  In  this  section  we  will  consider  the  case  in  which 
e»d.  and  p are  real.  The  case  of  complex  constants  is  considered  in' the  next 
section. 

First,  we  may  set  o'.  - l because  it  is  merely  a scale  factor. 
Second,  the  magnitude  of  (3  must  not  be  greater  than  unity  for  W^(s)  is 
then  unstable.  If  |p|  ^ 1,  the  typical  term  is  stable. 

To  determine  the  loci  of  typical  terms  of  the  form  (3-ll),  we  need 
only  apply  some  of  the  rules  of  the  loci  of  complex  functions.  First,  note 
that  the  locus  of  1 ♦ {3  e-*^  is  a circle  of  radius  |p|  centered  at  the 
point  (1.0).  To  find  the  locus  of  W^(jco),  the  inverse  of  a circle  must  be 
found.  If  Ifl I *1,  the  circle  (l  ♦ e"^)  passes  through  the  origin, 
and  its  inverse  is  a straight  line  parallel  to  the  imaginary  axis.  If 
IpI^CI,  the  circle  does  not  pass  through  the  origin,  and  its  inverse  is 
another  circle.  The  loci  of  two  stable  first-degree  partial  fractions  are 
shown  in  Fig.  3»2. 
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3.23  Second-Degree  Partial  Fractions  (complex  roots) 

If  a typical  term, 

V.(s) 2L- , (3-12) 

1 1 ♦ p e 

•sT 

in  the  partial  fraction  expansion  of  a rational  function  of  e has  complex 
constants,  there  will  be  another  term, 

rJ/# 

W2ts)  - — » (3-13) 


in  the  expansion  whose  constants  are  the  conjugates  of  those  of  W^(s). 

(The  asterisk  means  conjugate.)  Both  W^(s)  and  W^Cs)  will  be  stable  if  and 
only  if  jpl^-1.  If  the  rational  function  has  real  coefficients  (as  in 
practical  problems)  the  terms  such  as  W^s)  and  W^Cs)  must  come  in  pairs. 
Add  the  two  and  obtain  a typical  second  degree  partial  fraction  with  real 


coefficients . 

W3(s)  - W^s)  ♦ W2(s) 


U ♦cxi*)  ♦ (uf* 

1 ♦ (P  ♦ P*)  e"8*  ♦ PP*  e'*®T 


(3-ll*) 


In  this  section  the  discussion  is  restricted  to  second-degree  partial  fractions 

-sT 

whose  denominators  have  complex  roots  of  e • To  simplify  the  analysis, 

W3(s)  can  be  written  as, 

1 ♦ 

W«(s)  - 

3 1 ♦ bxe 


V 


-sT 


-sT 


♦ b.e 


-2  at 


(3-15) 


which  differs  from  (3— lU)  by  only  a constant  multiplying  factor. 

To  insure  that  the  roots  of  the  denominator  of  (3-15)  are  complex, 
o 

b^  4 Ubg*  With  this  condition  imposed  on  the  constants  of  the  denominator 
of  (3-15),  it  can  be  considered  a basic  building  block  of  a program  transfer 


function 
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Now  let  us  consider  the  stability  of  such  a building  block. 

W^(s)  will  be  stable  if  and  only  if  both  W^(s)  and  Wg(s)  are  stable.  A 
comparison  of  (3-lii)  and  (3-15)  shows  that, 

bx  - p ♦ p*,  and  b2  - PP*  - |p|2.  (3-16) 

Therefore,  the  necessary  and  sufficient  condition  for  the  stability  of 
W^(s)  is  that,  0^  b^  ^1.  We  may  combine  this  with  the  condition  for 
complex  roots  in  the  denominator  of  W^(s)  and  obtain, 

(rj  < b?  < 1 

as  the  necessary  and  sufficient  condition  that  insures  stability  of  W^(s) 
and  complex  roots  of  its  denominator. 

In  Fig.  3.3  there  are  plotted  three  loci  of  building  blocks 
of  the  form  (3-15).  All  three  are  stable.  It  should  be  observed  that  by 
changing  only  the  numerator  of  the  transfer  function,  three  completely 
different  loci  have  been  obtained. 

By  adding  building  blocks  of  the  form  discussed  in  these  section*;, 
a desired  frequency  characteristic  can  be  synthesized.  The  synthesis  of 
a differentiating  program  is  discussed  in  a subsequent  section. 
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3.3  Realization  of  Programs  From  Their  Transfer  Functions 

The  purpose  of  this  section  is  to  develop  and  compare  methods  by 

-aT 

which  a program  can  be  derived  from  a rational  function  z * e which  is 
the  transfer  function  of  the  program.  How  this  rational  function  is  arrived 
at  in  the  first  place  is  the  concern  of  the  last  section  of  this  chapter. 

3.31  General  Considerations  in  Program  Realization 

In  choosing  a particular  method  of  programming,  one  may  consider 
the  following  factors!  storage  requirements  and  time  requirements.  To  a 
certain  extent  one  of  these  requirements  can  be  reduced  at  the  expense 
of  increasing  the  other  and  the  optimum  method  will  depend  on  the  particular 
application.  It  is  necessary,  therefore,  to  make  available  various  possible 
methods  of  programming  and  to  form  some  idea  about  the  requirements  of  each; 
intelligent  program  realization  can  then  be  adapted  to  each  application. 

In  the  consideration  of  storage  requirements  of  linear  programs, 
it  is  convenient  to  distinguish  three  types  of  storage*  data  storage,  constant 
storage,  and  instruction  storage^".  The  data  are  the  successive  sampled 
values  of  input  and  output.  The  complexity  of  a program  is  closely  related 
to  the  number  of  constants  and  to  the  age  of  the  data  to  which  the  program 
refers.  The  progran  can  be  divided  into  arithmetic  and  manipulative  parts. 

The  number  of  arithmetic  operations  involved  is  roughly  proportional  to 
the  number  of  constants,  each  implying  a multiplcation  (of  a piece  of  data 
by  the  constant)  and  an  addition  (of  the  product  to  the  other  terms). 

The  number  of  manipulative  operations  is  related  to  the  "age"  of  the  oldest 

I 

It  is  understood  that  in  a general-purpose  computer  there  is  no  physical 
difference  between  the  storage  registers  containing  numbers  or  instructions, 
and  any  register  may  hold  either  kind  of  information.  The  distinction 
made  here  is  only  for  the  purpose  of  discussion. 
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data  used,  where  age  is  expressed  in  terms  of  sampling  intervals.  All 
"younger"  data  must  also  be  stored  even  if  not  used  at  each  calculation 
(the  corresponding  constants  being  zero),  for  eventually  they  will  become 
the  oldest  data.  After  each  calculation  of  a new  output  value,  the  manipulative 
instructions  shift  each  piece  of  data  to  a storage  location  at  which  an  older 
piece  of  data  has  been,  the  oldest  data  being  lost.  The  manipulative  program 
is  seen  to  rearrange  the  data  storage  in  such  a manner  that  at  the  new  sampling 
point  the  same  arithmetic  program  will  calculate  a new  output  value . 

The  time  requirement  of  a program  i3  the  product  of  the  number  of 
instructions  to  be  carried  out  and  the  average  duration  of  an  instruction. 

The  latter  factor  depends  on  the  physical  characteristics  of  a particular 
computer  and  is  more  or  less  fixedj  the  number  of  instructions  performed 
in  sequence,  however,  depends  in  part  on  the  manner  of  program  realization. 

In  each  particular  realization  a significant  trading  of  time  for  storage 
is  possible  by  so-called  cyclic  procedures.  One  notes  that  often  the  calculation 
of  each  term  in  a program  involves  the  same  sequence  of  arithmetic  operations. 

The  simplest  and  fastest  procedure  is  to  store  as  many  of  these  sequences 
as  there  are  terms  to  be  calculated.  Considerable  storage  may  be  saved, 
however,  by  storing  these  instructions  only  once  and  cycling  through  them  as 
many  times  as  there  are  terms  to  calculate.  Unfortunately,  the  time  requirement 
increases  considerably,  for  in  each  cycle  the  addresses  of  the  instructions 
must  be  adjusted  to  make  them  refer  to  different  storage  locations  for  different 
terms  and  the  number  of  cycles  must  be  counted  to  permit  termination  of 
the  cycling  process. 

The  following  sections  and  related  appendices  will  serve  as  specific 
illustrations  of  these  considerations  in  programming. 
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302  Direct  Regression  or  Direct  Programming 

The  starting  point  of  our  realization  procedure  is  the  general 


expression  for  the  transfer  function  of  a linear  program, 

-msT 

j 

■=nsT  * 


-sT  -2aT  , -msT 

a ♦ »e  ♦ a.e  ♦ ...  ♦ a e 
01  2 n 

W(s)  - ^ 


1 ♦ be  8 ♦ b2e 


(3-18) 


♦ . . . ♦ b e 
n 


- ...  - b^(t  - nT), 


In  order  to  interpret  a program  in  the  time  domain  it  is  necessary  to  eliminate 

fractional  expressions.  The  most  straightforward  way  of  doing  this  follows 

directly  from  (3-18).  From  it  we  can  obtain 

tf(a)  - (a-  ♦ ...  ♦ a e""8*)  r(s)  - (to. e”aT  ♦ ...  ♦ b e"MT)  tf(s). 
um  x ' n 

(3-19) 

The  inverse  transform  of  this  expression  is 

to(t)  - a^i(t)  ♦ a^itt  -?)♦...♦  a^Ct  " ®T)  - to^o(t  - T) 

(3-20) 

where  o(t)  and'i(t)  are  impulse-modulated  (sampled)  time  functions  having 

Ibhe  value  zero  everywhere  except  at  the  sampling  points.  In  terms  of  some 

continuous  functions  o'(t)  and  i(t),  which  agree  with  the  area-values  of 

cf(t)  and  ^(t)  at  the  sampling  points,  (3-20)  is  often  written  as 

o.  * a i.  ♦ a_i.  ,♦...  + ai.  - b.o.  . - .;.  - bo.  , (3-21) 

j o j l j-1  m j-m  1,  j-1  n j-n* 

where  j signifies  particular  sampling  point  and  j-k  the  k-th  preceding  sampling 

point.  Eq.  (3-21)  is  more  familiar  to  the  numerical  analyst  than  (3-20), 

but  the  two  are  entirely  equivalent  and  are  called  regression  formulas. 

These  equations  state  that  the  present  result  (output)  is  conputed  by  a 

finite  linear  combination  of  the  present  and  past  input  values  and  of 

past  results  (output  values).  * 

Several  characteristics  of  regression  formulas  should  be  observed. 

If  the  right  side  of  (U-20)  or  (i*-2l)  has  at  least  one  non-zero  b^,  then 


1 
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the  present  output  depends  on  at  least  one  previous  output,  which  in  turn 
depends  on  an  output  further  back  and  so  on.  It  follows  that  the  present 
output  value  is  affected  by  output  values  as  far  back  as  the  start  of  the 
problem  and  therefore,  also  by  input  values  that  far  back.  Thus  regressing 
to  a finite  number  of  output  values  corresponds  to  regressing  to  an  unlimited 
number  of  input  values.  This  aspect  of  the  regression  equations  is  important 
and  will  be  further  emphasized  in  the  following  sections. 

Interesting  conclusions  can  be  drawn  concerning  memory  requirements 
on  the  digital  computer  by  considering  the  actual  programming  of  the  regression 
formula,  (3-21).  Assuming  that  none  of  the  coefficients  a^.  and  b^  are  zero, 
one  can  easily  see  that  in  order  to  calculate  a new  output  value  o^,  when 
a new  input  value  i^  is  received,  m previous  input  values  and  n previous 
output  values  will  have  to  have  been  remembered,  requiring  m ♦ n memory 
positions  for  data  where  m and  n are  the  subscripts  of  the  last  non-zero 
coefficients,  and  furthermore,  it  is  necessary  to  store  all  these  data  even 
if  some  other  coefficients  are  zero  because  at  the  next  sampling  point 
the  same  pieces  of  data  will  be  associated  with  different  coefficients. 

Pt  can  be  stated  that,  at  least  when  programming  is  done  by  the 
illustrated  direct  regression  method,  the  data  memory  consists  of  m ♦ n 

registers  (memory  positions)  where  m and  n are  the  degrees  of  the  numerator 

~sT 

and  denominator  polynomials  in  z ■.  e of  the  program  transfer  function  W(s). 
Actually  this  data  storage  requirement  may  be  reduced,  as  will  be  shown  in 
Sections  3.33  and  3*3k> 

To  be  able  to  make  comparisons  between  the  various  synthesis 
procedures  it  is  necessary  to  do  the  actual  programming.  This  exercise 

i 

is  left  to  the  appendices,  and  the  results  will  be  compared  after  the 
other  synthesis  procedures  will  have  been  discussed.  In  Appendix  B the 
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arithmetic  and  manipulative  parts  of  the  direct  regression  program  are  first 
constructed  separately;  then  a new  more  compact  program  is  shown  which  inter- 
leaves the  arithmetic  and  manipulative  instructions.  Although  the  Whirlwind 
code  is  used,  the  results  and  conclusions  can  be  considered  quite  general 
in  view  of  the  fact  that  the  instruction  complements  of  most  general-purpose 
digital  computers  are  conspicuously  similar* 

3*33  Cascade  Programming 

-sT 

If  the  numerator  and  denominator  polynomials  in  z ■ e are 
factored,  (U-18)  takes  the  form 


W(s)  - a 


o 


(1  ♦ c^e”8^)  (l  ♦ Cge”87)  ...  (l  ♦ cme_8T) 

Cl  ♦ dje"8*)  (1  ♦ dge”87)  ...  (l  ♦ d^"8*)  * 


(3-22) 


where  -(l/c^)  and  -(l/d^)  are  the  roots  of  numerator  and  denominator 
respectively,  when  considered  as  polynomials  in  z*  Because  the  coefficients 
a^  and  b 

or  will  come  in  conjugate  pairs.  At  any  rate,  it  is  possible  to  group  the 
monic  factors  of  (3-22)  in  some  manner 

W(s)  - W^s)  W2(s)  ...  Wp(s),  (3-23) 

where  each  W^(s)  is  of  a rational  form  in  z having  a numerator  and  a 
denominator  of  not  higher  degree  in  z than  W(s)  itself  has. 

The  form  of  (3-23)  reminds  one  of  the  transfer  function  of  cascaded 
linear  units  in  a servo  system.  Cascading  means  that  the  output  of  one 
unit  becomes  the  input  to  the  next  one.  There  is  no  difficulty  in  using 
the  same  interpretation  to  define  cascaded  programs.  At  every  sanpling 
point  the  output  of  each  regression  equation  is  used  in  calculating  the 
output  of  the  next  one* 


k of  these  polynomials  are  real,  the  c^  and  d^  will  also  be  real 
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To  be  more  specific,  let  us  assume  that:  (l)  W(s)  is  a proper 

rational  fraction  in  z,  that  is,  m<  n^;  (2)  all  roots  -l/c^  and  -l/d^ 

are  real  and  distinct,  since  generalization  to  the  case  of  conjugate  complex 

roots  turns  out  to  be  direct:  (3)  a - 1 in  order  to  avoid  its  nuisance 

o 

2 

value  in  the  discussion  . With  these  assumptions  it  is  possible  to  have 
p * n in  (3-23)  with  the  denominator  of  each  W^(s)  being  a single  monic 


factor  in  if  that  is. 


- - OJ 

1 ♦ c,  e 

V> 

1 ♦ V 


(3-2U) 


where  may  or  may  not  be  zero,  but  d^  / 0.  There  are  n factors  of  the 
type  (3— 2U)  each  representing  a simple  regression  equation.  In  m of 
the  factors  c^  / 0,  in  the  other  n-m  factors  c^  ■ 0.  The  data  storage 
associated  with  each  Wk(s)  equals  2 when  c^  / 0,  and  1 when  c^  ■ 0.^ 
However,  the  input  data  that  must  be  stored  when  c^  / 0,  is  also  the  output 
data  that  had  to  be  stored  for  the  preceding  cascade  program  ^(s)j 
consequently,  there  is  only  one  data  to  be  stored  for  each  W^(s)  regardless 
of  the  value  of  c^,  except  for  the  first  one  W^(s) . But  when  m < n (proper 
rational  fraction),  one  ck  - 0,  say  c^  ■ 0,  making  the  total  required  data 


If  m^n,  W(s)  can  be  written  as  the  sum  of  a polynomial  and  a proper 
rational  fraction  in  z - e"s^.  The  program  corresponding  to  the  polynomial 
part  is  a simple  linear  combination  of  input  values.  Discussion  of  this 
case  is  omitted  without  any  serious  loss  of  generality. 

^ If  / 1,  only  a simple  multiplication  has  to  be  added  to  the  program. 

^ Each  W^(s)  is  the  transfer  function  of  just  a regression  equation  and 

its  data  storage  is  the  sum  of  degree  of  numerator  and  denominator,  as 
discussed  in  the  previous  article. 
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storage  n;  a material  reduction  over  the  m ♦ n data  needed  in  direct  regression 
programming. 

In  order  to  translate  the  cascade  scheme  into  an  actual  program, 
we  may  proceed  as  follows.  First,  we  write  (3-23)  (with  p * n)  in  terms 
of  input  and  output  transforms,  as 


%)  ^(b)  0-(s)  0 (s) 

____  ■ 1 e C • Ml  a H 

T(s)  ^(s)  T2( s)  T(s) 


(3-25) 


One  way  of  making  (3-25)  an  identity  is  by  letting 

\M  - ?( s) 

?2(s)  - Q^(s) 

I3(s)  - ^(f) 

■ Vi(s) 

which  make 

0(s)  - 'O  (s). 

n 

Using  the  relations  (3-25)  and  (3-26)  we  obtain 
Ots)  'O  (s)  'ol(s)  %) 

■ ± • c e mi  • ___ 

iM  Tf(s)  'ol(s)  '(T  . (s) 

i n— i 


(3-26) 


(3-27) 


The  rarious  factors  equal  the  respective  program  transfer  functions;  namely. 
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0^(8) 

W1(s)  - 

1 

•\ 

7(a) 

1 ♦ 

. -ST 

V 

v>  . 

W2(s)  - 

e 

1 ♦ 

c2.->T 

<^(a) 

1 ♦ 

d2e 

1 

V (3-28) 

0(.) 

• 

• 

v*>  • 

1 ♦ 

-sT 
c e 

n 

Vl<*> 

1 ♦ 

, -sT 
d e 

n 

— J 

wherfc  m of  the  c^  are  not  zero*  Multiplying  by  the  denominators  changes 
the  set  C3-28)  into 

(X  ♦ 4l«'“T)'o1(>)  -7(a) 

Cl  ♦ dj.-”1)  tf/a)  - (1  ♦ o2e-aT)  O^a) 

(1  ♦ dj.'*1)  03(a)  - (1  * c3a'sT)  V?(a) 

e 

U ♦ dne”ST)  " (1  ♦ cn*”aT)  ^V-I(s) 

The  inverse  transform  of  the  foregoing  set,  with  one  term  of  each  equation 
transposed  to  the  right  side,  is  the  desired  set  of  regression  equations* 
OjCt)  -T(t)  - d^o^(t-T) 

oTgCt)  - ol(t)  ♦ o2o^(t-T)  - djS^Ct-T) 
o^(t)  * o^(t)  ♦ c3o£  (t-T)  - d^Ct-T) 

e 

^(t)  - o'  . (t)  ♦ c ol  (1-T)  - dftt-T) 

n-x  n n-i  n 


(3-29) 


(3-30) 
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The  detailed  coded  program  corresponding  to  (3-30)  is  shown  in  Appendix  C. 

Cascade  programming,  although  not  referred  to  by  that  name,  is 
a familiar  technique  in  numerical  procedures.  However,  the  clear-cut  and 
general  equivalence  of  the  direct  regression  and  cascade  programming  is  not 
always  well  understood.  Cascade  programming  arises  naturally  from  the  kind 
of  thinking  prevelant  in  numerical  work.  Consider  the  simple  example  of 
solving  the  second-order  differential  equation. 


♦ Py  - 0. 

dt 


(3-31) 


The  derivatives  may  be  considered  as  the  separate  variables,  y'  (t)  and 
y"( t)j  then  we  obtain  the  following  three  sampled  functions* 
y"(t)  - -py(t  - T) 

T'(t)  • ^"(t)  ♦ y'(t  - T)  l (3 


'y(t)  - Ty’(t)  *y(t  - T) 


(3-32) 


where  the  first  equation  of  the  set  is  derived  from  (3-31)  while  the  second 
and  third  are  elementary  first-difference  extrapolations.  The  set  (3-32) 
indicates  cascade  programming  because  the  output  of  the  first  equation  is 
in  the  input  of  the  second,  and  the  output  of  the  second  equation  is  the 
input  to  the  third.  The  peculiar  thing  in  this  case  is  that  the  input 
to  the  first  equation  is  not  an  independent  function  but  directly  related 
to  the  output  of  the  last  equation.  This  feature  establishes  the  constraint 
imposed  by  the  differential  equation. 

The  Laplace  transform  of  the  set  (3-32)  is 

t*  - - e"sTr 

¥•  - TV*  * e"sTY  ( (3-33) 

Y - iT  ♦ e“aTf*  I 


(3-33) 
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from  which  the  explicit  relations  between  inputs  and  outputs  are  obtained 
as  follows  t 


T"  - - e "sTI 


T'  - 


1 - e 

_X 


-sT 


Y» 


1 - e 


For  realization  by  three  cascaded  factors,  we  have 


W1(s)  - -pe 

W,(s) X 

£ 1 - e 

W,(s) 31 


-sT 


-sT 


1 - e 


-sT 


(3-3U) 


(3-35) 


It  is  clear  that  a single  transfer  function  can  be  made  to  replace  the 
cascaded  system  of  three;  thus 

W(s)  - WL(s)W2(s)W3(s) 


V(s)  - 


-pT2e~sT 


(3-36) 


1 - 2e’8T  ♦ e‘2sT 


The  corresponding  regression  equation  is  simply  obtained  as 

T(s)  - (2  - pT2)e"sTY(s)  - e”2sTY(s).  (3-37) 

The  inverse  transform  of  this  equation  is 

y(t)  - (2  - pT2)y(t  - T)  -J(t  - 2T).  (3-38) 

which  could  have  been  otained  from  (3-32)  by  the  elimination  of  y' (t)  and 
yN(t),  but  even  in  this  simple  case  the  process  of  elimination  in  the  time 
domain  is  not  direct* 
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The  fact  is  that  in  numerical  work  a cascade  method  such  as  (342)  ie 
much  more  generally  used  than  the  direct  regression  of  (3-38)#  Often  there 
is  good  justification  for  this  preference;  for  instance  the  values  of  the 
first  and  second  derivatives  may  also  be  needed.  However,  when  such  or 
similar  justifications  do  not  exist,  the  direct  regression  may  turn  out  to 
be  simpler  than  cascading.  In  the  present  example,  (3-32)  calls  for  one 
more  constant,  two  more  multiplications  and  one  more  addition  than  (3-38). 

If  the  first  two  equations  of  the  set  (3-32)  are  combined,  one  multiplication 
is  saved;  furthermore,  the  manipulations  in  the  direct  method  happen  to  be 
more  awkward.  Because  in  this  case  the  input  and  output  are  the  same  quantity 
the  formulas  of  Appendices  B and  C are  not  directly  applicable,  the  requirements 
of  the  two  methods  must  be  determined  by  actual  trials. 

3.3U  Parallel  Programming 

If  the  transfer  function  of  a program  is  expanded  by  partial 
fractions  in  terms  of  z,  (3-18)  tales  the  form 


VCs)  « 


f 

n 

1 ♦ d e 
n 


(3-39) 


as  long  as  m n.  Thus,  the  traisfer  fraction  W(s)  is  replaced  by  the 
sum  of  a number  of  simpler  transfer  functions;  namely, 

W(s)  - v^Cs)  ♦ w2(s)  ♦ ...  ♦ wp(s),  (340) 

where  some  of  the  W.  Cs)  may  be  the  combination  of  several  partial  fractions, 
but  all  are  of  lower  degree  than  W(s)  itself. 

The  form  of  (340)  may  remind  one  of  parallel  combinations  of 
network  admittances.  Paralleling  means  that  the  same  input  (driving  voltage) 
is  applied  to  all  component  admittances  and  the  output  (driving-point  current) 
is  obtained  as  the  sum  of  individual  outputs  (current  through  each  admittance)* 
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The  same  interpretation  can  be  applied  to  parallel  programming.  The  programming 
mill  involve  p regression  equations  all  using  the  same  input  values,  and  all 

I 

their  outputs  adding  to  produce  the  desired  over-all  output. 

To  arrive  at  a more  specific  interpretation,  me  first  make  a few 
restrictions  again:  (l)  W(s)  is  a proper  fraction,  i.e.,  m < n,  and  (2) 

the  roots  of  the  denominator  polynomial  are  real  and  distinct.  Then  all 
constants  f^  and  d^  of  (3-39)  are  real  and  in  (3-liO)  p can  equal  nf  moreover, 
each  term  of  (3-39)  is  a simple  regression  equation  involving  two  constants 
and  one  data  storage.  Thus,  the  total  number  of  data  to  be  stored  is  only  n. 
Just  like  in  the  case  of  cascade  programming,  the  lower  requirement  for  data 
storage  of  parallel  programming,  may  be  a great  advantage  over  the  direct 
programming  method.  However,  this  feature  does  not  mean  that  parallel  or 
cascade  programming  should  always  be  enployed  in  preference  to  direct 
programming.  For  instance,  there  is  the  case  when  m ■ 0}  i.e.,  the  numerator 
of  W(s)  is  1 (or  aQ).  Of  the  input  values  the  program  uses  only  the  present 
one  and  the  total  data  storage  is  n regardless  of  the  programming  scheme 
usedj  on  the  other  hand,  the  number  of  constants  will  be  n for  the  direct 
and  cascade  method,  but  2n  for  the  parallel  method,  putting  the  latter  at 
a disadvantage.  Similarly,  if  the  denominator  of  the  over-all  transfer 
function  lacks  several  terms  (say,  the  denominator  is  1 - bne”nsT),  then 
factorization  of  the  denominator  introduces  all  terms,  making  the  cascade 
and  parallel  program  much  longer  than  the  direct  program.  Another  factor 
which  may  militate  against  the  use  of  parallel  programming  is  the  presence 
of  multiple  roots  in  the  denominator.  If  a root  is  of  multiplicity  r, 
it  may  produce  up  to  r terms  of  degrees  r,  (r-l),  ...2,1  (the  r-degree 
term  never  being  absent)  in  the  partial  fraction  exapnsion,  but  the  same 
root  will  require  only  one  r-degree,  or  r first-degree,  cascaded  factors. 
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In  order  to  interpret  the  parallel  method  of  programming,  we 
proceed  in  the  usual  manner.  For  the  various  terms  of  (3-1*0 ) with  p ■ n, 
we  write 

_ £■» 


and 


wnC.) 


W(s) 


'l(a) 

1 ♦ d^8* 

V 

f(s) 

• 

i ♦ v"*T 

• 

0(e) 

n 

f 

, n 

I(s)  1 ♦ d e 
n 


ff(s) 


-=bT 


J 


1(a) 


(3-ia) 


(3-1*2) 


Cross-multiplication  by  the  denominators  in  (3-1*1)  yields  the  set 
(1  ♦ dLe'8T)01(s)  - f^s) 

(1  ♦ d2e'8T)02(s)  -f2I(s) 


(3-1*3) 


(1  ♦ d e‘lT)0  (s)  - f ill)  J 

n n n — ^ 

while  in  view  of  (3-1*1)  and  (3-1*2),  (3-1*0)  can  be  written  as 
ote)  - 0^(s)  ♦ 0^(s)  ♦ ♦ (T (s)  • 


(3-1*1*) 
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The  inverse  transforms  of  (3-^3)  and  (3-Ui)  yield  the  desired  set  of 
regression  equations,  which  follows. 


3^(t)  - f^Ct)  - dL^t-T) 

^2(t)  - f2i(t)  - dgtfgCt-T) 

o' (t)  - ft(t)  - dSl(t-T) 
n n n n 

3(t)  -^(t)  ♦3'2(t)  ♦ ...  ♦ ^n(t) 


(3-U5) 


The  detailed  coded  program  corresponding  to  (3-U5)  is  shown  in  Appendix  D. 

Parallel  programming  has  not  been  generally  used  in  numerical 
work.  To  the  knowledge  of  the  writer,  the  usual  methods  of  numerical 
analysis  do  not  naturally  lead  from  a direct  regression  equation,  which 
has  reference  to  several  previous  input  and  output  values,  to  a set  of 
simpler  regression  equations,  each  of  which  refers  only  to  the  last 
input  value  and  to  a preceding^-  output  value.  By  the  method  of  frequency 
transformation  the  parallel  method  is  found  quite  directly. 


In  case  of  complex  d^'s  in  (3-39),  a combination  of  two  conjugate  complex 

partial  fractions  in  z will  result  in  a slightly  more  complicated  regression 
equation,  involving  one  additional  input  and  output. 
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3.35  Comparison  of  Programming  Methods 

The  purpose  of  this  section  is  to  compare  the  effectiveness  of 
the  various  methods  of  program  realizations  based  on  the  transfer 
functions  of  the  programs.  A complete  general  treatment  appears  too 
far-fetched  and,  therefore,  this  study  is  limited  to  a certain  class  of 
programs.  Despite  these  limitations,  which  are  discussed  below,'  the 
investigation  is  sufficiently  general  to  show  how  the  results  can  be 
used  to  improve  the  instruction  code  of  a general-purpose  computer  or 
to  design  a special-purpose  computer,  when  these  are  used  in  control 
applications. 

The  three  methods  which  will  be  compared  are  listed  belowt 

(a)  direct  programming, 

(b)  cascade  programming, 

(c)  parallel  programming. 

Other  programming  schemes  may  be  derived  from  the  rational  transfer 
function  W(s).  One  may. carry  out  the  long  division  in  z of  the  numerator 
by  the  denominator  until  he  arrives  at  a certain  number  of  terms  of  the 
quotient.  The  transfer  function  can  then  be  expressed  as  the  sum  of  the 
quotient  terms  and  of  the  remainder  divided  by  the  divisor  (the  original 
denominator) . Any  number  of  variations  can  be  obtained  by  stopping  the 
long  division  after  different  number  of  steps,  but  only  in  the  most  unusual 
cases  can  this  approach  be  expected  to  yield  a more  efficient  scheme  of 
programming  than  the  three  major  methods  discussed  in  the  preceding  sections. 
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Other  schemes  that  are  even  more  artificial  than  the  long-division  scheme  may 
be  derived,  but  no  other  general  programming  method  has  been  found  that  gives 
promise  of  effectiveness  comparable  to  the  three  which  are  considered.  It 
is  noted  that  in  certain  cases  a combination  of  two  of  the  three  listed 
methods  may  turn  out  to  be  more  efficient  than  any  one.  An  example  of 
such  a case  is  described  below. 

As  the  basis  of  comparison  of  programming  methods,  the  requirements 
in  storage  and  time  are  used.  The  particular  application  or  purpose  decides 
which  of  these  two  factors  should  deserve  more;  attention.  It  is  assumed 
that  the  complete  sequence  of  instructions,  as  used  at  each  sampling  point, 
is  stored;  the  possibility  of  cycling  programs,  which  re-uses  a short 
sequence  of  instructions  for  the  calculation  of  each  term,  is  not  discussed. 
Essentially  the  Whirlwind  I code  is  used  throughout,  but  variations  are 
considered. 


As  a starting  point  we  recall  that  the  transfer  function  of  a 
linear  program  is, 

-sT  -2sT  . -msT 

* *V  *V  * — * V (3-1*6) 


W(a)  - 2^- 
r(s) 


...  *V-MT 

1°.  lyr-sT  » b2.-Sst  * ...  * b 


n 


In  general,  m and  n may  be  any  positive  integeps^  and  indeed,  their  relative 
sizes  will  hardly  influence  the  comparisons  to  follow.  Nevertheless,  it  is 
helpful  to  distinguish  three  cases: 

~sT 

(l)  n » 0.  (3-U6)  reduces  to  a polynomial  in  z - e~  ; i.e., 

the  new  output  value  depends  only  on  present  and 
past  input  values,  not  on  past  outputs  also. 
Present-day  numerical  analysis  abounds  in  numerical 


t 


This  is  in  contrast  with  networks  where  certain  restrictions  on  the 
degress  of  numerator  and  denominator  polynomials  often  exist. 


■f 
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•*> 


processes  corresponding  to  this  special  case* 

(2)  ■ < n*.  (3-46)  has  the  form  of  a proper  rational  function 

of  z in  this  case*  In  Sections  3*32,  3*33,  and  3*34 
dealing  with  the  various  programing  schemes,  this 
case  was  assumed  for  the  sake  of  simplicity. 

(3)  it  ^n.  The  rational  function  in  z of  (3-46)  maybe  called 

improper,  but  it  can  be  converted  to  the  sum  of  a 
polynomial  (Base  l)  and  of  a proper  rational  fraction 
(Case  2)  in  z - e . 

In  order  that  the  storage  and  time  estimates  to  be  arrived  at  should  apply 
to  all  cases,  it  is  necessary  to  define  the  following  quantities  with 
reference  to  (3-46)* 

m * degree  of  numerator, 

- n - degree  of  denominator. 


"k 


one  less  than  the  number  of  non-zero  constants  in 
the  numerator  (ra^^  n). 

one  less  than  the  number  of  non-zero  constants  in 
the  denominator  (n^  ^ n) . 


m ■ one  more  than  the  excess  of  m over  nj  i.e.. 


m ■ m-n+1  when  m is  n,  and  m - 0 otherwise.  For 
e ^ e 

proper  rational  fractions  m < n and  m ■ OC 
On  basis  of  the  coded  programs  shown  in  the  appendices,  the  table  of 
Fig.  3*4  summarizes  the  storage  and  time  requirements  in  terms  of  the 
quantities  just  defined.  This  tabulation  is  more  general  than  the  results 
given  in  the  appendices,  for  in  the  appendices  it  was  also  assumed  that 
none  of  the  constants  were  zero,  that  is,  m^  * m and  n^  - n;  furthermore. 


Examples  are  numerical  methods  based  on  polynomial  approximations  with 
equidistant  spacing  of  the  independent  variable.  Indeed,  such  examples 
form  not  an  insignificant  portion  of  the  available  numerical  techniques* 
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only  case  (2)  was  treated  making  re  - 0.  On  the  other  hand,  in  the  tabulation 
of  Fig.  3.k  these  restrictions  of  the  appendices  are  absent,  but  the  following 
assumptions  are  still  made:  the  roots  of  the  numerator  and  denominator 

are  real  and  distinct,  and  the  straightforward  programming  techniques  of 
the  appendices  is  used.  Thus,  the  constants  stored  are  those  that  appear 
explicitly  in  the  various  regression  equations.  Actually  some  saving  in 
instructions  would  result  from  the  use  of  certain  ratios  of  these  constants. 
For  instance,  the  regression  equatibn  [cf.  (B-U 5)J 

^(t)  - f-jitt)  - d^Ct-T)  (3-U7) 

takes  six  instructions,  as  shown  in  the  coded  program  of  Appendix  D. 

If,  however,  (3— U7)  is  written  as 

Sj(t)  - fx  |j(t)  - , (3-U8) 

its  coding  would  cost  five  instructions  only,  but  certain  questions  on  the 
relative  sizes  of  the  constants  would  arise.  It  seemed  best  to  avoid  such 
questions,  because  the  considerations  here  are  rather  general  and  the  value 
of  a too-specialized  treatment  is  questionable. 

The  comparison  of  the  three  methods  of  programming  can  be  undertaken 
by  considering  each  item  of  Fig.  3»U»  Because  of  the  straight  sequential 
programming  the  time  requirements  are  the  same  as  the  storage  for  instructions 
and,  therefore,  consideration  of  storage  will  give  a complete  picture. 

As  far  as  the  number  of  constants  stored  are  concerned,  the  direct 
method  is  not  worse  than  the  cascade,  which  in  turn  is  not  worse  than  the 
parallel  method.  This  is  so  because  in  the  direct  method  only  the  non-zero 
constants  of  (3-U6)  have  to  be  stored,  while  factorization  in  the  cascade 
case  will  produce  as  many  constants  as  there  are  roots  in  z.  In  the  parallel 
method  two  constants  (root  and  residue)  are  produced  for  each  denominator 
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( ) 


( ) 


root  and  if  the  numerator  is  not  of  lesser  degree  than  the  denominator, 
further  terms  and  constants  result.  As  an  example,  consider 


V(s) 


1 

F 


-lisT 

6 


C3-U9) 


for  which 

m - 1,  - 1,  \ m 0 

n - k,  \ m 

According  to  the  table  of  Fig.  3«U  the  various  constant  storage  requirements 
are 


direct!  m^  ♦ n^  ♦ 1 » U 

cascade:  m ♦ n ♦ 1 » 6 

parallel:  2n  * 8 

These  figures  can  be  simply  checked.  In  the  direct  case  the  four  constants 
are  apparent  in  (3 *4j9).  For  the  cascade  case,  the  transfer  function  is 
written  as 

5 _ 1 e-*T 

w(s)  • x — ” — T — IgT > 

l*ie  11  1-i*'1  1 ♦ — e sr  1 - — e ar 

n ft  ,2  2 


and  the  six  constants  in  question  are:  ♦l/^T,  -1  /[Tt  +l/2,  -l/2,  *5/8, 

and  -l/h.  The  manner  of  programming  illustrated  in  Appendix  C actually 
necessitates  the  separate  storing  of  positive  and  negative  constants, 
even  though  6f  the  same  magnitude* 

For  parallel  programming  W(s)  of  (3-U9)  is  expanded  in  partial 
fractions  in  terms  of  z and  takes  the  form 
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V(b) 


5 ♦ 2l/2~ 


x'fT 


r^r 


5 - 2/r~~ 
8~ ^ 

rrr^" 

/r 


1 ♦ ^ e 


9 

IS" 

i~-  8r  * 


i 

’IT 

: — i — =5T“ 
1 - 2a 


(3-51) 


The  eight  constants  to  be  stored  are  evident  in  the  foregoing  filiation. 

The  next  item  of  comparison  is  the  data  storage,  for  which  the 
above  example  reads,  on  basis  of  Fig*  3.U 
direct*  m ♦ n - 5 
cascade*  n * k 

parallel*  n ■ U 

These  figures  can  be  verified  in  the  three  foregoing  equations.  The  numerator 

of  (3-U9)  indicates  that  one  past  input  value  (corresponding  to  the  e term) 

must  be  stored;  the  present  input  is  used  as  it  arrives  and  then  stored 

as  the  past  input  for  the  next  calculation,  as  shown  in  Appendix  B.  Thus 

the  numerator  implies  one  data  register  only.  Similarly  the  denominator 

-sT 

implies  the  storage  of  four  past  output  values,  even  though  the  e and 
-3sT 

e terms  are  absent;  for  the  corresponding  past  outputs  must  be  remembered 
for  the  next  calculation. 

For  the  cascade  method,  (3-5>0)  seems  to  indicate  five  past  data 

-sT 

to  be  remenbered;  however,  the  e term  of  the  last  numerator  refits  to 
a past  input  that  is  also  the  past  output  of  the  preceding  factor,  since 
in  cascade  programming  the  input  of  a component  program  is  the  output  of 
the  previous  one* 

In  case  of  parallel  programs  the  four  past  data  are  quickly 
identified  with  the  e terms  of  (3-51)* 


1 
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The  expressions  for  instruction  storage  and  for  time  requirement* 
are  identical,  and  produce  the  following  tally  in  the  present  exanqplet 
directs  2(m  ♦ m^  ♦ n ♦ n^)  ♦ 7 ■ 23 
cascades  3»  ♦ Un  ♦ 6 ■ 25 

parallels  7n  ♦ 1*  - 32 

No  verification  of  these  figures  is  carried  out  by  detailed  coding  of  the 
programs  because  the  appendices  cover  the  general  case.  The  advantage 
seems  to  be  on  the  side  of  direct  programming  as  far  as  time  is  concerned, 
but  this  advantage  is  slight  and  arises  from  the  fact  that  in  the  present 
example  two  denominator  constants  are  zero.  An  advantage  of  direct  programming 
appears  also  in  the  total  storage  requirements  for  the  same  reasons 
directs  1*  ♦ 5 ♦ 23  -32 

cascades  6 + l*  + 25  ♦l*  36 

parallels  8 + l*  + 32  + 1*  1*5 

This  example,  as  well  as  the  tabulation  of  Fig.  3 .1*,  indicates 

the  disadvantage  of  parallel  programming.  It  seems  that  this  kind  of 
programming  may  have  an  advantage  over  either  of  the  other  two  in  certain 
cases,  but  hardly  ever  over  both  at  the  same  time.  Thus,  the  choice  narrows 
down  to  direct  and  cascade  programs,  or  possible  combinations  thereof. 

To  show  how  a combination  of  methods  may  be  used,  we  write  (3-1*9)  as 


V(s) 


1 

l~=2sT 


5 1 -sT 

ff  “IT* 

l~-TiT 


(3-52) 


1 ’ f 6 1 " g e 

which  indicates  a cascade  combination  of  two  direct  programs,  for  which 
respectively 


*1* 
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■ - 0 

v 0 

n - 2 
\‘1 


m ■ 1 

“k  " 1 
n - 2 

nk  - 1 


m ■ 0 ra  * 0 

e e 

The  direct  program  of  each  cascaded  component  is  somewhat  simpler  than  it 
would  be  for  two  separate  direct  programs  because  the  input  and  output  devices 
are  manipulated  only  once  for  the  composite  program,  rather  than  once  for 
each  component  program.  This  saving  amounts  to  six  instructions,  tbaa  the 
instruction  storage  or  time  requirement  is: 

first  components!  2(m  ♦ m^  ♦ n ♦ n^)  +6  - 12 

second  components!  2(m  ♦ m^  ♦ n ♦ n^)  ♦ 6 * 16 
saving  as  indicated  above  -6 

* 

total  instructions  '22 


Four  constants  appear  in  (3-*»2),  two  of  which  are  accidentally  identical, 
and  one  of  which  is  made  1;  thus,  the  constant  storage  isi 


first  component: 

- 2 

saving 

- 

1 

second  component! 

- 3 

saving 

- -1 

2 

total  constants  3 

A saving  arises  in  data  memory  also,  because  the  past  input  of  the  second 
component  is  also  the  past  output  of  the  first  one.  This  gives  the  following 
need  of  data  storage  > 
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first  component:  m ♦ n ■ 3 

second  conponent:  m + n • 3 

saving  ^1 

U 

The  results  of  this  example  are  summarized  in  Fig.  3.5,  which  shows  a small 
advantage  of  the  mixed  method  over  the  direct  one. 

To  pursue  further  the  detailed  comparison  of  these  various  methods 
of  programming  would  lead  to  undue  specializations  in  the  Whirlwind  code  and 
to  results  of  doubtful  general  value.  The  illustrated  attack  on  the 
realization  problem,  however,  shows  how  a useful  estimate  of  the  complexity 
of  coded  programs  can  be  gained  from  the  evident  properties  of  their  transfer 
functions.  Three  further  problems  will  be  touched  on  briefly:  (l)  computing 

delays,  (2)  means  of  using  the  results  to  select  computer  codes;  and  (3) 

t 

means  of  using  the  results  to  design  special-purpose  computers. 

A consideration  that  has  been  omitted  in  our  discussion  is  the 
delay  incurred  through  the  computation  itself.  If  a digital  computer  is 
used  as  part  of  a number  of  control  systems  — say,  50  systems  — , then 
in  each  sampling  interval  it  performs  50  computations,  one  for  each  system. 

The  time  of  a computation  is  then  at  most  l/$0  of  the  sampling  time,  T, 
and  this  delay  is  presumably  negligible.  If,  however,  one  digital  computer 
were  used  with  each  system,  the  computation  may  and,  for  the  sake  of  efficiency, 
should  take  an  appreciable  part  of  the  sampling  time.  Such  a delay  would 
be  very  serious  and  the  computer  would  have  to  perform  a prediction  in 

addition  to  the  required  compensation.  In  turn,  this  would  lengthen 

/ 

the  program,  make  it  less  effective,  and  may  even  force  a longer  sampling 
time;  indeed,  in  a marginal  case,  in  which  the  original  compensating  program 
had  a delay  nearly  as  large  as  the  sampling  time,  the  effect  may  become 
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cumulative,  since  a longer  sampling  time  would  in  turn  require  a better 
and  longer  program,  and  so  on.  In  such  marginal  cases  and  in  any  case  in 
which  the  computing  time  is  not  negligible  with  respect  to  the  sampling 
time,  the  direct  programming  has  a tremendous  advantage  over  all  other 
methods.  A glance  at  the  direct  regression  equation  (3-20)  shows  clearly 
that  ail  terms  but  the  first  one  on  the  right  side  of  the  equation  can 
be  computed  before  the  new  input  value  is  obtained.'*'  The  computing  delay 
will  thus  be  the  time  of  merely  calculating  the  term,  aQT(t),  and  adding 
it  to  the  already  prepared  partial  result.  This  delay  may  conceivably 
be  negligible* 

All  realizations  of  real-time  linear  programs  involve  accumulation 
of  products  as  their  arithmetic  action  and  the  transfer  of  data  from  one 

2 

register  to  another  as  their  manipulative  action.  In  case  of  a single-address 
instruction  code,  such  as  that  of  Whirlwind,  the  ex  (exchange)  operation 
was  shown  in  Appendix  B to  be  very  helpful  in  improving  the  efficiency 
of  the  code.  Other  Improvements  axe  possible  by  incorporating  special 
operations  which  facilitate  the  particular  type  of  programs  on  hand. 

Computers  using  multiple-address  codes  could  be  particularly  efficient  in 
such  applications.  For  instance,  in  a three-address  code  an  instruction 
could  locate  a constant,  a piece  of  data,  and  transfer  that  data  to  a 
third  address,  after  which  it  would  multiply  the  constant  and  data 

^ The  second  composite  program  In  Appendix  B is  written  in  this  manner* 

2 

Each  instruction  specifies  an  operation  and  the  storage  address  of  a 
single  operand. 

3 

This  operation  exchanges  the  contents  of  the  accumulator  register  with 
the  specified  storage  register.  Thus,  one  instruction  performs  a double 
duty  by  obtaining  hew  data  from  storage  and  also  transferring  to  storage 
a partial  result* 
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accumulating  this  product  with  the  partial  result  always  left  in  the  arithmetic 
element  of  the  computer.  This  single  order  would  complete  both  the  arithmetic 
action  (accumulation  of  products)  and  the  manipulative  action  (transfer  of 
data  to  an  "older-data"  register)  associated  with  one  terra  of  a regression 
equation. 

Similar  considerations  allow  one  to  adapt  special-purpoee  or 

fixed-program  digital  computers  to  control  specifications.  To  be  somewhat 

specific  we  assume  that  the  computer  is  used  as  part  of  a single  control 

system  and  will  have  to  perform  only  one  computation  in  each  sanpllng  period. 

The  computer  would  not  operate  appreciably  faster  than  one  computation  per 

sampling  period  and  in  order  to  minimize  the  computational  delay  it  would 

follow  a direct  regression  program.  In  order  to  keep  such  a single-system 

computing  equipment  from  becoming  excessive,  a serial1  computer  would 

probably  be  used.  The  program  of  the  computer  would  be  fixed  to  correspond 

to  a direct  regression  program  of  certain  complexity  as  defined  by  the 

degrees  of  m of  the  numerator  and  n of  the  denominator  of  the  program 

transfer  function.  The  constants,  could  be  set  manually  on  toggle  switches 

2 

or  relays,  or  they  could  be  stored  on  the  same  high-speed  storage  device 
on  which  the  data  are  stored.  A serial  adding  unit  with  proper  switching 
equipment  would  allow  the  multiplication  of  constant  and  data  (by  repeated 
additions)  and  the  addition  of  such  product  to  the  accumulated  partial 
result.  The  physical  size  of  such  a digital  control  unit  maybe  quite 
feasible  in  certain  applications  and  the  design  of  such  a simple  special-purpose 

digital  computer  would  be  particularly  Justified  if  the  incoming  data  were 

» 

sampled  and  digital  to  start  with. 

A serial  computer  operates  on  each  digit  of  a number  in  sequence;  thtoa, 
the  equipment  is  not  duplicated  for  each  digit. 

2 

Magnetic-drum  memory,  for  instance. 
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3«U  Synthesis  of  Programs  In  the  Frequency  Domain 

3.1il  General  Synthesis  Procedure 

The  synthesis  of  computer  programs  in  the  frequency  domain  may 

be  broken  down  into  the  three  following  stages  (l)  specification  of  the 

desired  frequency  characteristic  or  lous  of  W(je&),  (2)  approximation  of 

the  desired  locus  by  a realizable  program  transfer  function,  and  (3) 

realization  of  the  program.  One  way  to  determine  the  desired  locus  is 

from  the  Laplace  transform  of  the  operation  the  computer  is  to  perform. 

The  second  step  is  the  difficult  part  of  the  problem.  The  desired  frequency 

-sT 

characteristic  must  be  approximated  by  a rational  function  of  e .No 
general  rules  are  available  for  making  this  approximation,  but  before 
making  the  approximation,  one  should  gain  some  experience  in  analyzing 
program  building  blocks  in  the  complex  plane.  Possibly  the  most  systematic 
approach,  at  present,  to  the  approximation  problem  is  to  make  successive 
approximations  to  the  desired  characteristic,  using  the  basic  program 
building  blocks  of  Section  3.2.  The  third  step  involves  only  a straight- 
forward inverse  Laplace  transform.  As  an  example  of  program  syhthesis 
in  the  frequency  domain,  a program  for  differentiation  will  now  be  synthesized. 
3.U2  Synthesis  of  a Differentiation  Program 

An  ideal  differentiator  establishes  the  following  relation  between 
input  and  output: 

o(t)  - ^ i(t).  (3-52) 

Disregarding  initial  conditions,  the  Laplace  transform  of  (3*52)  is 

H(.)  - 2W  . (3.53) 

«») 

and  this  is  the  desired  transfer  function.  For  s ■ jco,  H(s)  becomes 
H(j«)  ■ J«. 


C3-5U) 
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So  the  locus  of  the  desired  transfer  function  is  the  imaginary  axis*  This 
completes  the  first  step  of  the  synthesis  procadure. 

-sT 

The  second  step  is  to  find  a rational  function  of  e that 
approximates  this  locus.  This  approximation  is  to  be  made  by  geometric 
considerations  based  on  the  desired  locus.  In  this  particular  example  it 
is  also  possible  to  employ  analytic  considerations  based  on  the  desired 
transfer  function  of  (3-5U)*  It  so  happens  that  in  the  present  case  the 
analytic  approach  is  simpler;  nevertheless,  the  geometric  approach  is  shown 
first. 

The  crudest  numerical  approximation  to  a first  derivative 
i$>  the  first  divided  difference. 

'S'(t)  ■ Jll-1).  (3.55) 

T 

The  Laplace  transform  of  (3-55)  is 

ots)  -T(.)  1 ~ »**  . (3-56) 

T. 

Thus  the  transfer  function  of  the  differencing  process  is 

w («)  ■ - 1 y~*T  . (3-57) 

Fig.  3*6  shows  the  locus  of  Wq( jco)  and  compares  it  with  the  desired  one* 


Figure  3*6  Comparison  of  first  derivative  and  first  difference  operators 
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with  respect  to  the  sampling  frequency)  the  two  loci  agree  reasonably  well* 


If  we  could  straighten  out  the  circular  locus,  we  would  have  a better 
approximation  of  the  desired  locus.  Figure  3*7  illustrates  a geometric 
construction  that  straightens  out  the  locus  and  gives  us  ideal  phase 


characteristics . 

A 21-  e"sT 


di(t-T) 

dt 


1 - e 


1 ♦ e 


(a)  Frequency  domain 


i(t)  - l(t-T) 
T 

di(t) 

- dt 


(b)  Time  domain 


Figure  3*7  Derivation  of  an  ideal  phased  realisable  differentiation 
operator 

The  vectors  (l/T)(l  - e-^00^)  and  (l/2)(l  ♦ e”*^)  are  drawn  for 
a particular  frequency*  Using  the  geometric  rule  that  a triangle  inscribed 
in  a semicircle  is  a right  triangle,  one  can  readily  show  that  M * W 
add  up  to  90°.  However,  since  o<  is  a positive  phase  angle,  it  must  be 
subtracted  from  0 (which  is  negative)  to  give  a resulting  angle  of  -90°, 
which  is  the  phase  of  an  ideal  integrator.  It  follows  that  division  of 
the  P -locus  by  the  o( -locus  will  yield  an  ideal-phase  formula*  The 


resulting  transfer  function  is 


W2(e) 


2 1 - e 

f : 

x 1 ♦ e 


-=Sf  9 


(3-58) 
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which  has  the  desired  phase  in  the  range,  0 4La>T  The  interpretation 

in  the  tine  domain  is  both  plausible  and  illunlnating.  The  inverse  transfom 
of  (3-58)  shows  that. 


ret)  ♦ y(t  - T)  _ i(t)  -T(t-T) 
2 


(3-59) 

(3-59)  states  that  the  average  of  the  derivatives  at  two  neighboring 

points  is  approximately  equal  to  the  divided  difference  for  those  points* 

It  is  interesting  to  note  that  the  sane  approximate  transfer 

function,  (3-58),  can  be  obtained  analytically  baaed  on  a rather  good 

-sT 

approximation  for  e • 


_aT~I  - | sT 
s — 

i ♦ |>t 


(3-*>) 


Solving  (3-60)  for  s yields 


8 ’ T 


2 1 - e 


-sT 


1 ♦ e 


(3-61) 


Although  in  this  particular  case  the  above  analytic  approach  is 
simple  and  fairly  accurate,  its  general  use  has  certain  drawbacks.  The 

most  obvious  one  is  that  the  rational  function  of  "s"  to  be  approximated, 

/ 

which  in  the  present  case  is  "s”  itself,  is  in  many  cases  not  explicitly 

known*;  rather  it  may  be  obtained  as  an  approximation  to  a desired  locus 

or  amplitude  and  phase  response.  Then  to  approximate  the  rational 

function  of  *s",  which  itself  is  but  an  approximation,  by  a rational 
-sT 

function  of  e puts  the  designer  on  shaky  grounds,  and  it  might  lead  to 
far  more  involved  programs  than  necessary.  There  is  no  substitute  to 
going  back  to  the  original  specifications  and  designing  directly  on  their 
basis.  Another  disadvantage  of  the  above  analytic  approach  is  that  it 
is  not  general.  One  could  replace  all  "s"  by  the  approximation  (3-6l), 
but  how  one  would  get  a better  solution  is  not  obvious* 
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We  have  an  approximation  of  the  differentiation  operator,  so  the 


next  thing  to  do  is  see  how  good  it  is.  Since  the  desired  locus  and  its 
approximation  lie  along  the  same  path,  a locus  study  does  not  give  a good 
comparison.  In  such  a case  separate  amplitude  and  phase  plots  can  be  studied. 
For  s ■ jco,  WjjCs)  becomes 


„ tjQ)  2 1 - .->T 


T . J2  . oaT 

3T  “ 'T  ttn  T » 


(3-62) 


which  verifies  the  previous  statement  that  W-(jco)  has  ideal  phase  characteristics. 


Hence,  it  is  sufficent  to  study  the  amplitude  characteristic  only. 

H(jco)  - jay 


therefore. 


w20) 


co T 

tan~2~ 


2 


(3-63) 


(3-6U) 


Thus  we  see  from  (3-61*)  that  the  ratio  of  the  approximate  function  to  the  ideal 
one  Is  always  greater  than  unity.  Figure  3.8,  a plot  of  the  amplitude 
characteristics,  shows  us  that  for  low  values  of  ©T,  say  for  <»T  » the 


| (1)  HC>) 

| (2)  W2(» 

| (3)  W-O) 


- 2 tan 

» 1.81]*  tan 


son  of  amplitude  characteristics  of  Differentiat 


rators 
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differentiation  program  will  give  quite  good  accruacy.  For  certain  control 

V If 

applications,  values  of  a>T  up  to  j or  even  ^ might  give  acceptable  accuracies. 

An  examination  of  the  amplitude  characteristic  of  W^Cs)  in  Fig.  3*8 
reveals  that  if  W^(s)  is  multiplied  by  a constant,  which  is  slightly  less  than 
unit,  we  will  obtain  a better  derivative  on  the  average.  The  new  transfer 
function  is 


W3(s)  - C W2(s)  - 

Let  us  arbitrarily  choose  C so  that 
2 C tan  ^ ■ j ; 


2 1 - e 

f 

1 ♦ e 


-sT 

=3T  * 


W^Cs)  ■ H(s)  for  sT 


(3-6$) 


J j.  Then, 

(3-66) 


so 


, 6 tan  £ 

The  improved  transfer  function  is 


0.907. 


(3-67) 


W3(s) 


1.8XU  1 - e 


-sT 


(3-68) 


1 ♦ e 

and  its  amplitude  characteristic  is  also  shown  in  Fig.  3.8. 

Both  curves  2 and  3 of  Fig.  3.8  accentuate  high  frequencies  which 
may  be  present  at  the  input  because  of  noise.  In  this  case,  a transfer  function 
whose  amplitude  characteristic  is  like  that  of  curve  h would  be  a more  desirable 
approximation  for  differentiation. 

The  inverse  transform  of  (3-68)  completes  the  synthesis  of  a 
differentiation  program.  The  result  is 


ott)  - [T(t)  - i(t-T)J  - o(t-T). 


(3-69) 


The  accuracy  of  this  differentiation  program  may  be  determined  from  Fig.  3*8, 
Curve  3 • 
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CHAPTER  IV 

FREQUENCY  ANALYSIS  OF  SOME  NUMERICAL  INTEGRATION  FORMULAS 

In  this  chapter  we  apply  the  methods  of  frequency  analysis  to 
several  numerical  integration  formulas:  the  trapezoidal,  Simpson's  1/3 

rule,  Simpson's  3/8  nils,  and  Ifeddle's  rule.  Frequency  analysis  is 
applied  to  determine  the  stability  of  these  formulas,  compare  their  ac- 
curacy, and  compare  their  transfer  functions  with  that  of  the  ideal 
integrator. 

U.l  Numerical  Integration 

In  the  numerical  integration  of  definite  integrals,  the  range 
of  integration  is  divided  into  a convenient  number  of  equal  intervals, 
and  the  values  of  the  integrand  are  defined  only  at  the  ends  of  these  inter- 
vals. Essentially  this  is  the  same  as  sampling  (or  impulse  modulating) 
the  integrand.  Let  the  distance  between  samples  be  T.  To  obtain  an 
approximate  value  of  the  integral  mb  may  determine  an  n^1  order  polynomial 
that  passes  through  n + 1 of  the  sampled  points  and  integrate  the  poly- 
nomial over  the  corresponding  range,  repeating  the  process  until  the  com- 
plete range  of  the  original  integral  has  been  covered.  If  the  sampled 
points  are  joined  by  straight  lines,  (approximation  by  a first  order  of 
polynomial)  the  re  silting  integration  formula  is  know  as  the  trapezoidal 
rule  (each  interval  of  the  integrand  is  approximated  by  a trapezoid). 

Joining  the  points  in  each  group  of  three  sampled  points  by  a parabola 
leads  to  an  integration  formula  known  as  Simpson's  l/3  rule.  If  the 
points  in  each  group  of  four  sampled  points  are  joined  by  a cubic  curve,  we 


Heport  B“225f 


-71- 


get  Simpson' s 3/8  rule.  The  trapezoidal  rule  and  Simpson's  1/3  rule 
are  quite  widely  known  and  used,  but  there  is  another  one,  called  Weddle's 
rule,  that  is  used  to  obtain  great  accuracy.  Joining  the  points  in  each  group 
of  seven  sampled  points  by  a sixth  order  polynomial  leads  to  Weddle's  rule. 

In  each  case  the  range  of  integration  should  be  divided  into  an  integral 
multiple  of  ”n"  intervals.  For  example,  to  use  Weddle's  rule,  the  range 

of  integration  should  be  divided  into  6,  12,  18 equal  intervals. 

In  what  follows  we  shall  designate  the  transfer  function  of  an 
ideal  integrator  as  H(s)s  Thus, 

H(s)  = i (U-l) 

with  which  the  approximate  integration  formulas  will  be  compared, 
lull  Trapezoidal  Rule 


Using  the  trapezoidal  rule  the  definite  integral. 


o(t) 

may  be  approximated  by. 


L 


i(x)  dx  • 


5l(t)  * | /[  i(t)  ♦ i(t  - T)  ] + [i(t  - T)  + i(t  - 2T) J + 

^ Cii-3) 


The  Laplace  transform  of  ( 1*— 3)  is 


0,(s) 


l 


Cl  ♦ .-*r)  (i  + .-*T)  ♦ .-2sr  ♦ 


] «»). 


(WO 
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.?> 


So  the  transfer  function  is 

°l(»)  m 

Wl(s)  “ TfsT  = 7 


1+  e 


-sT 


1 - a 


-=Sf 


(k-$) 


A little  algebra  shows  that  for  s = 


■TOT  " r 


»T  obT 

COt  9T" 


(U-6) 


U.12  Simpsons  1/3  Rule 

Using  Simpson’s  l/3  rule  the  definite  integral  (U-2)  may  be 
approximated  by 


[t 


°2Ct)  3 5 1 Li(t>  + ' T> 


+ i(t  - 2T)]  + 


[i(t  - 2T)  + U(t  - 3T)  + i(t  - UT)  ]+ 


(V?) 


The  Laplace  transform  of  (li-8)  is 


o2Cs)  = * i(.)  (i  + U e*sT  + . -2,TJ  (l  + .-2,T  ♦ . -UT  +, 

CI.-6) 


or. 


_ T 


®2^S^  3 5 


, ^ ) -sT  , -2sT 
1 + 4 e + e 


1 - e 


-=7ST 


I(s) . 


(k-9) 
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Therefore,  the  transfer  function  for  Simp  a on' s l/3  rule  is 


W2(s) 


Vs)  _ T 1 + It  e“sT  + e -2,T 


TT«T  I 


1 - e 


“=2sT 


(4-10) 


Dividing  (4-10)  by  (4-1),  letting  s = J®  and  using  some  algebraic  and 
trigonometric  manipulations  leads  to  the  ratio 


W2^)  _ »T  2 ♦ cos  rnT 

BT3®7  T sin  »f 


C4-n) 


( ) 


4*13  Simpson's  3/8  Rule 

The  approximation  to  the  definite  integral  (2*— 2)  that  is  obtained 
using  Simpson's  3/8  rule  is 

o3(t)  * ^i(t)  «■  3i(t  - I)  ♦ 3i(t  - 2T)  + i(t  - 3T)  J ♦ 


E 


i(t  - 3T)  ♦ 3iCt  - UT)  ♦ 3i(t  - $T)  ♦ i(t  - 6t) 


The  Laplace  transform  of  o^(t)  is 


) 

(4-12) 


,(.)  * * !(•)  (l  + 3 . -T  ♦ 3 .-2*1  ♦ *'3‘T  + •"6ST+ 


or. 


0,(s) 


3T  1 + 3 • x + 3 • * » 


-2sT  A -3sT 


1 - a 


■^TT 


2 I(  s) 


(4-13) 


(4-14) 


f > 


1 
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Hence  for  Sin?) son's  3/8  rule,  the  transfer  function  is 


W,(s) 


03(e) 

TTiT 


3T  1 + 3 e’sT  + 3 e“2aT  ♦ e“3BT 


1 - e 


(U-l$) 


For  s » jca,  the  ratio  of  to  H(j»)  is 


. 3.T 

1®  - T" 


1 + cos  coT 


(1+2  cos  ottT)  tan  ^ 


(U-16) 


A considerable  amount  of  manipulation  is  required  to  obtain  the  above  form. 


U o lU  Weddle r s Rule 


E|y  Weddle's  rule  the  approximation  of  the  definite  integral 


(lt-2)  is 


«!.(*) 


s{D 


i(t)  + 5i(t  - T)  + i(t  - 2T)  + 6i(t  - 3T)  + 


i(t  - UT)  + £i(t  - $T)  + i(t  - 6T)J  + |i(t  - 6T)  + 

5i(t  - 7T)  + i(t  - 8T)  ♦ 6i(t  - 9T)  + i(t  - 10T)  + 

5i(t  - 11T)  + i(t  - 12T)  1+ I 

J j (U— 17) 

In  the  same  manner  as  before,  the  transform  of  o^(t)  is 


y*> 


3T  1 ♦ 3 ♦ .~2,T  ♦ 6 .'3,T  ♦ .-1*1  ♦ S.**  ♦ .-6,T 


1 - e 


I(s); 


CU— 18) 


Report  R-225 


-75- 


eo  the  transfer  function  for  Weddle's  rule  is 


V") = 


V“> 

ttst 


3T  1+5  .-*T  + .'2,T  ♦ 6.*3,t  + .-U,T  + S.-5aT+  .-6”1 

to  

(U-19) 


By  using  a considerable  amount  of  algebraic  and  trigonometric  manipulation, 
we  get  for  s = 


V^Cjea)  3 1 3 cos  a>T  + cos2  «T 

0(3<a7  *"  5 (1+2  cos  »T)  sin  o»T  * 


(U-20) 


Uo2  Comparison  of  Numerical  Integration  Formulas 

With  the  above  equations,  we  can  get  a complete  picture  of  the 
four  approximation  formulas  in  both  the  time  and  frequency  domains. 

Equations  (k-5),  (U*  10),  (U-15),  and  (U-19)  are  the  transfer  functions 
of  each  of  the  numerical  integration  processes  and  from  these  the  stability 
of  each  one  can  be  determined.  Let  us  now  examine  the  denominator  of  each 
transfer  function.  If  the  change  of  variable,  s * e , is  made,  it  is 
easily  seen  that  the  magnitude  of  the  roots  of  all  the  denominator  poly* 
nomials  is  unity}  however,  there  are  no  multiple  roots.  Therefore,  each 
of  the  numerical  processes  is  stable. 

Now  we  must  consider  the  accuracy  of  each  of  the  integration 
formulas.  Equations  (U-6),  (U-ll),  (U-16),  and  (U-20)  give  the  ratio  of 
the  particular  transfer  function  to  that  of  the  ideal  integrator.  In  Figure 
U.l  these  ratios  are  plotted  as  functions  of  o&T,  and  we  see  clearly  that, 
of  the  four,  Simpson's  1/3  rule  and  Weddle's  rule  are  the  best  for  a>T 
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below  -■  radians.  For  example,  suppose  that  we  wish  to  integrate  a sine 
wave  of  radian  frequency  and  want  the  error  to  be  less  than 
For  each  of  the  approximation  formulas,  how  many  samples  mast  be  taken 
in  a cycle  of  the  sine  wave?  The  answer  can  be  obtained  rapidly  from 
Figure  li0l  by  noting  the  frequencies  at  which  the  amplitude  ratios  beoome 
0.975  or  1.02$,  as  listed  below. 


Trapesoidal 
Simpson's  1/3  Rule 
Simpson's  3/8  Rule 
Weddle's  Rule 


«0T  = 30°  ; 12  sample s/cycle 
»0T  * 75°  ; 5.8  " ■ 

® T = 6o°  y 6.0  » » 

O 

CO  T = 80°  J u.5  " " 

o 


The  number  of  samples  per  cycle  is  indicated  for  each  rule,  and  thiB  is 
obtained  by  dividing  360° 


by  the  indicated  angle. 
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APPENDIX  A 


Proof  that  the  Locus  of  QCjaa)  Croasea  the  Real  Axis  either  Normally 
or  Tangentially  at  to  ■ 0 and  HF 


Recall  that  Q(s)  is  given  by  the  polynomial  in  e 


-at 


ja.  . ...  _ 

QCa)  - > b.  e"8T  ; b - 1 . 
k - 0 * 0 ■ 

For  e-O  and  — QCa)  ia  real  becauae  each  term  of  the  polynomial  ia 

real.  Since  the  locus  ia  symmetrical  about  the  real  axis,  it  must  cross  the 
real  axis  at  these  points* 

In  order  to  examine  the  behavior  of  the  locus  of  Q(ja)  at  these 
point 8,  take  the  derivative  of  Q(s)  with  respect  to  a* 


k - 1 


e-T 


dQ 


-at. 


Observer  that  ^ is  also  a polynomial  in  e f therefore,  it  will  also  be 
real  for  s - b and  — J ^ . 


Now  consider  the  derivative  in  the  neighborhood  of  s - o and 
♦ J . If  / 0 at  these  points,  we  will  prove  that  the  Q-locus  crosses 

the  real  axis  perpendicularly.  In  the  region  of  interest  let  ds  - J $ , 

where  S is  a small  increment  of  oo.  Since  must  be  real  (and  unequal 

to  zero  as  we  have  assumed),  dQ  * )dQ  j in  order  to  this  condition. 


(Q.E.D.) 

d©  4 >, 

We  must  now  discuss  the  case  in  which  ■gj  ■ o for  a - 0 or  -Jajr. 
dQ  1 

First1  observe  that  if  ^ - 0,  Q must  have  a saddle  point  in  the  region 

dQ 

near  the  point  where  ^ - 0.  Let  us  now  make  the  change  of  variable, 

I 7 — r 

For  an  excellent  discussion  6f  the  behavioul*  of  functions  near  saddle 
points,  see  Quillemin,  "The  Mathematics  of  Circuit 
and  Sons,  New  Tork,  19U9,  pp.  298-302* 
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«=>sT 

z * e , so  that  Q(s)  becomes  D(z)  which  is 

m 

D(b)  - 5 


TT=T> 


bk  ! • 


In  the  immediate  vicinity  of  a saddle  point,  the  function  behaves  as 

D(z)  C ♦ C (z  - z )P 
op  o 

in  which  the  C’s  are  constants,  z is  the  value  of  z at  which  the  saddle 

* o 

point  occurs,  and  "p  - 1"  is  the  order  of  the  saddle  point.  In  this  case, 
zq  * +1.  In  plotting  the  locus  of  D(z),  we  map  the  unit  circle  of  the  z-plane 
into  the  D-plane  (see  Fig.  A-l). 


Consider  the  map  in  the  vicinity 
of  a possible  saddle  point  (z  ■ *l). 
Observe  that  for  z near  z 

O' 

z - zq  - dz  ♦j/dzj.  There- 
fore, in  the  vicinity  of  a saddle 

point  D(z)  is 

D(z)  - Cq  ♦ Cp  (♦;)  |dz  |)P. 


Fig  o A-l  Unit  circle  in  the  z-plane  that 
maps  into  the  D-plane 

This  readily  shows  that  if  "p"  is  even,  the  locus  in  the  D-plane  (or  Q-plane) 
is  tangent  to  the  real  axLs.  If  "p"  is  odd,  the  locus  is  normal  to  the  real 
axis. 

We  will  now  summarize  the  results  obtained. 

a)  If  / 0 for  s ■ 0 or  the  locus  of  Q(joj)  is  normal  to 

the  real  axis  for  s » 0 for  ■♦j ^ respectively. 
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b)  If  ~ - 0 for  a - 0 or  ♦ Q(s)  has  a saddle  point  at  the 

-sT 

point  where  the  derivative  is  zero.  The  change  of  variable,  z ■ e , 

permits  us  to  write,  D(z)  - C ♦ C (z  - z for  z in  the  immediate  vicinity 

o p o 

of  the  saddle  point.  If  "p"  is  even,  the  Q-plane  locus  is  tangent  to  the 
real  axis  at  the  saddle  point.  If  "p"  is  odd,  the  locus  is  normal  to  the 


real  axis 
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APPENDIX  B 

Coding  of  Direct  Regression  Programs 


The  regression  formula 

'fc(t)  - a T(t)  ♦ a-T(t-T)  ♦ ...  ♦ aT(t-mT)  - b.'Stt-T)  - ...  - b'oTt  -nt) 
ox.  ra  X n 

(B- 

is  to  be  programmed.  Assume  that  the  data  and  the  constants  are  stored 
as  follows t 


Register 


Content 

(Constants) 


Register 

No. 


Content 

(Data) 

!<t) 

i(t  - 

• 

T) 

• 

• 

TCt  - 

mT) 

?r(t  - 

T) 

^t  - 
• 

2T) 

• 

• 

'oCt  - 

nT) 

Partial  and 
final  results 

These  registers  are  not  used  in  the  second  composite  program 
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The  program  will  first  be  coded  in  two  distinct  parts*  arithmetic  and 
manipulative.  The  arithmetic  part  performs,  at  each  sampling  point,  the 
arithmetic  operations  called  for  by  the  above  regression  formula  and  thus 
calculates  a new  output  value* 


Register 

No. 


First  Program,  Arithmetic  Portion J 


Content 

(Instruction) 

Result 

ca  O.n 
mr  B.n 

ts  R.O  

V -b  ?F(t  - nT) 

^ n 

ca  O.n-1 ' 
mr  B.n-1 
ad  R.O 
ts  R.O 

I -b  U(t  - nT) 
n 

P.(hn-h) 

etc. 

1 

ca  0.1 
mr  B.l 
ad  R.O 

ts  R.O  — 

P.CW3) 

ca  I.m 
mr  A.m 
ad  R.O 

ts  R.O  

bk3tt  - kT) 


The  code  is  explained  in  Sc.D.  Thesis  "Treatment  of  Digital  Control  Systems 
and  Numerical  Processes  in  the  Frequency  Domain,"  J.M.  Salser,  Appendix  l.G 
Vpl.  2,  August  1,  1951*  M.I.T# 
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Continued* 


Register 

No. 

Content 

(Instruction) 

Result 

• 

• 

e 

e 

• 

0 

P.U(n*n) 

ca  1.0 

etc. 

mr  A.O 

ad  R.O 

ts  R.0  — 

-y  o(t) 

si 

selects  the  relevant  output  device 

(as  specified  by  the  address  section) 

P.(W*W5) 

rc  R.O 

records  output,  o(t),  into  output 
device 

It  is  clear  that  in  this  illustration  each  term  of  the  regression  equation 
costs  U instructions. 

After  the  calculation  of  6(t)  at  a particular  sampling  point 
the  data  storage  has  to  be  rearranged  for  the  next  calculation*  the  present 
■3{t  - nT)  can  be  lost,  all  other  “^Ct  - kT)  are  to  be  stepped  down  one  storage 
register,  and  the  new  output  value,  o'(t),  just  computed  is  put  into  0.1; 
the  rearrangement  of  the  iT(t  - kT)  is  analogous,  and  the  new  input  value  to 
be  received  goes  into  1.0.  The  coded  program  performing  these  manipulations 
follows. 


4 


P.(6*+6n+6) 


sp  P.l 


1 


-8U- 

HOGRAM,  MANIPULATIVE  PORTION 


Description 

/ moves  75 

J of 'o(t  - 

t - (n-l)T 
nT)  and.  Ic 

into  location 
>ses  o(t  - nT) 

moves  o 
J "3  (t  - (r 

t - (n-2)T 
~l)l]  local 

into 

Aon 

moves  b(t  - T)  into  U(t  - 2T) 
location 


moves  3(t)  into  o(t  - T) 
location 


moves  i jt  - (m-l)^  into  i(t-mT) 
location  and  loses  T(t-mT) 


moves  i(t)  into  i(t  - T)  location 


selects  the  relevant  input  device 
(as  specified  by  the  address  section) 
and  makes  computer  wait  until  device 
receives  a new  input  value 

reads  content  of  input  device  into 
i(t)  location 


returns  control  to  beginning  of 
whole  program. 
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The  manipulations  are  seen  to  cost  2 instructions  per  terra  of  the  regression 
equation. 

There  are  various  ways  in  which  this  program  can  be  streamlined. 
The  main  considerations  are  storage  and  time.  It  is  possible  to  save 
substantial  storage  (with  a sufficiently  long  regression  equation)  by 
programming  the  6 instructions  (U  arithmetic  and  2 manipulative)  required 
for  each  term  only  once  and  using  them  over  and  over  for  the  various  terms, 
each  time.  In  order  to  do  so,  a short  program  must  be  added  to  change  the 
appropriate  address  sections  in  the  6 instructions,  which  can  thus  be 
made  to  compute  a different  term  each  time.  This  address -changing  routine 
materially  lengthens  the  time  of  calculation,  unless  some  very  specialized 
instructions  or  equipment  is  designed. 

It  appears  more  desirable  to  concentrate  on  reducing  the  time 
requirements  in  most  control  applications.,  for  storage  is  easier  to  increase 
than  speed,  which  seems  to  be  the  ultimate  limitation  in  the  applicability 
of  digital  computers  to  controlling.  In  our  present  example  a notable 
reduction  in  time,  and  also  in  storage,  results  from  mixing  the  arithmetic 
and  manipulative  steps  and  using  a new  instruction,^  ex,  which  exchanges 
the  content  of  the  storage  register  specified  by  its  address  with  the 
content  of  the  accumulator.  The  corresponding  coded  program,  which  still 
uses  the  same  constant  and  data  storage,  follows. 


T 


This  instruction  is  actually  used  in  Whirlwind  applications  on  a temporary 
basis.  The  code  used  for  this  instruction  Is  £e  to  indicate  its  temporary 
nature;  final  adoption  of  this  instruction,  however,  is  likely. 
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( I 


Register 
' No. 


P.(l*n*l*m-8) 


Content 

(Instructions) 

♦ 


mr  A. 2 
ad  1.2 
ex  r.i 

ts  1.2 

nr  A.l 


ad  1.1 


ta  0.1 


si 


rd  1*1 

mr  A.O 
ad  0.1 
ts  0.1 

si  

rc  0.1 


P.(Un+Um+6) 


sp  P.l 


Result 


T 


into  Storage : partial  result 

into  AC i i(t  - T) 

T(t  - T)  into  T(t  - 2T)  location 


into  Storage!  partial  result 
(note  content  of  0.1  has  already- 
been  used  so  that  this  register 
is  available) 


selects  the  relevant  input  device 
(as  specified  by  the  address  section) 
and  makes  computer  wait  until  device 
receives  a new  input  value 

£e ads  content  of  input  device, 
i(t)  into  Itt  - T)  location 


into  Storage:  final  result 

tT(t)  into  7S(t  - T)  location 


selects  the  relevant  output  device 
(as  specified  by  the  address  section) 

records  output,  “ott),  into  output 
device 


returns  control  to  beginning  of 
program 
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The  above  composite  program  is  seen  to  result  in  considerable 

saving  of  storage  and  time  over  the  first  program,  which  was  given  mainly 

for  illustrative  purposes.  It  uses  four  instructions  per  terra  calculated 

rather  than  6,  and  even  saves  two  data  registers,  1.0  and  R.O.  Register  1.0 

is  not  needed  because  the  incoming  data  is  immediately  used  in  the  calculation 

while  register  R.O  is  superfluous  because  the  partial  results  can  be  stored 

in  the  register  from  which  the  data  has  just  been  removed  for  calculation. 

One  should  note  another  important  advantage  of  the  second  program! 

to  all  practical  extent,  it  eliminates  computational  delays  entirely. 

This  is  so,  because  all  the  computation  is  performed  in  advance  of  receiving 

the  input,  and  when  the  input  value  T(t)  is  received,  there  are  only  a few 

instructions  to  be  carried  out  in  order  to  obtain  the  output,  'STt).  Only 

direct  regression  programming  has  this  advantage. 

The  tally  of  direct  regression  composite  programming  in  terms  of 

m and  n,  the  degree  of  numerator  and  denominator  polynomials  of  the  program 

transfer  function,  is  as  follows: 

Time  requirement  (in  number  of 
instructions  to  be  carried  out 
in  sequence  at  each  sampling)  lm  ♦ lin  ♦ 6 


Storage  Requirements: 

Constants 

Data 

Instruction 

Total 


ra  ♦ n ♦ 1 
m ♦ n 

Urn  * lin  » 6 
6a  ♦ 6n  ♦ 7 


The  above  tally  is  made  under  the  assumption  that  none  of  the  constants 
are  zero.  If  some  constants  are  zero,  the  constant  and  program  storage, 
as  well  as  the  time,  requirements  will  be  reduced,  but  not  the  data  storage 
requirement.  These  more  specific  requirements  are  taken  into  account  in 
the  summary  of  Art.  3*35* 
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APPENDIX  C 


Coding  of  Cascaded  Programs 


The  set  of  regression  equations 
o^t)  ■ i(t)  -d^Ct-T) 

o2(t)  - o^(t)  ♦ c^Ct-T)  -d2o2(t-T) 

4 cn3n-l(t-T)-dn'i5(t-:I2l 


Cc-i) 


is  to  be  programmed.  Assume  the  following  arrangements  of  number  storage) 


egister 


ontent 
(Constants) 


ontent 
(Data) 


•^(t-T) 

^(t-T) 


o(t  - T) 


Partial  Result 


A.O 
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In  the  coded  program  to  follow  it  is  assumed  that  none  of  the  indicated 
c^  is  zero;  i.e.,  m - n - 1.  Variations  are  easily  accounted  for.  The 
program  instructions  follow. 


Begister 

No. 


Content 

(Instruction) 


rd  R.O 


ca  0.1 


rar  D.l 


ad  R.O 


ex  0.1 


mr  C.2 


ad  0.1 


ts  R.O 


ca  0.2 


mr  D.2 


ad  R.O 


Description 


selects  input  device  and  waits  until 
device  has  new  input  valuer  ±(t) 

reads  T(t)  into  temporary  location 


-dj3^(t  - T)  obtained 
TJ^(t)  obtained 


to  Storage*  1T-(t)  into 3-(t-T) 
location 

to  ACi  o^Ct  - T) 


o^(t)  ♦ “ T)  obtained 

to  Storage!  partial  result 


d^Ct)  obtained 
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Continued i 


Register 

No. 


F.(7n-S) 

etc. 


Content 

(Instruction) 


ex  O.n-1 


mr  O.n 
ad  O.n-1 
ts  R.O 
ca  O.n 


Description 

to  Storage:  7Tn_^(t)  into 

Z . (t-T)  location 
n-1 

to  AC:  Z - (t-T) 
n-1 

?5  , (t)  ♦ c Z , (t-l)  obtained 
n-1  n n-1 

to  Storage:  partial  result 


mr  D.n 
ad  R.O 
mr  A.O 


o(t)  obtained 


ts 


si 


O.n 


to  Storage:  c>(t)  into  Z( t-T) 

location 

select  output  device 


rc  O.n 


P.(7n*3) 


sp  F.l 


records  output,  Z( t),  into  output 
device 

returns  control  to  beginning  of 
program 


Thus,  if  n ■ n - 1,  the  program  is  7n  ♦ 2 instructions  long.  Suppose 
m ■ n - 2 and  let  c^  ■ 0;  then  the  sequence  P.6  through  P.12  above  would 
be  replaced  by  the  following  shorter  sequence  P.'  .9  through  P'  .12 . 


P‘  .9 

ts  0.1  — 

->  puts  fl^(t)  into  Tf^(t-T)  location 

P'  .10 

ca  0.2 

P»  .11 

mr  D.2 

P»  .12 


ad  0.1 


(S^Ct)  obtained 
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Rl 


Thus,  each  c^  ■ 0 saves  3 instructions. 

The  tally  for  cascade  programming  can  now  be  written: 
Time  Requirements:  3m  ♦ Un  ♦ $ 


m ♦ n ♦ 1 
n 

1 

3m  .♦  lm  ♦ 6 
iim  ♦ 6n  ♦ 8 


Comparison  of  these  requirements  with  those  of  other  methods  of  programming 
is  done  in  Art.  3*3!>» 


selects  input  device  and  waits 
until  device  has  new  input  value! 
i(t) 

reads  £(t)  into  its  assigned 
storage  register 

to  Storage:  f^itt) 

— to  AC i Z (t  - T) 

•»  to  Storage:  ^ (t) . 

-y  to  Storage:  f?T(t) 

— to  AC:  S^Ct  - T) 


to  Storage:  o^(t) 


to  Storage:  fjT(t) 

to  AC:  7>  (t  - T) 
n 
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Register 

No. 

Content 

(Instructions) 

Description 

mr  D.n 

a 

ad  O.n 

ts  O.n  — 

->■  to  Storage!  "o' (t) 
* n 

P.(6n+3) 

ad  O.n-1 

oft)  ♦ o'  . (t)  obtained 
n n-i 

etc. 

• 

ad  0.n*2 
• 

etc. 

• 

• 

ad  0.2 

ad  0.1 

7$(t)  obtained 

P.(7n+2) 

ts  1,0  — 

to  Storage*  75(t) 

etc. 

si  ___ 

selects  output  device 

rc  1.0 

records  output,  *8(t)  into 
output  device 

P.(7n*5) 

sp  P.l 

returns  control  to  beginning  of 
program 

In  parallel  programing  none  of  the  indicated  constanta  can  be 


zero,  and  the  only  possible  saving  is  when  several  constants  have  the  same 
value.  Even  then  the  program  itself  is  not  affected  materially. 

The  tally  for  parallel  programming  follows! 


Time  Requirement  7n  ♦ 5 

Storage  Requirements! 

Constants  2n 

Data  " n ♦ 1 

Instructions  7n  ♦ 5 

Total  lOn  + 6 


Further  discussion  of  these  requirements  is  left  to  Art.  3*35* 


